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CHAPTER  1: 
Introduction 


1.1  Motivation 

"Unleash  us  from  the  tether  of  fuel." 

-Lieutenant  General  Jamis  Mattis, 

Future  Fuels,  Naval  Research 
Advisory  Committee  Report 
April  2006 

During  the  past  decade  of  conflict,  service  members  on  bases  across  Iraq  and  Afghanistan 
could  hear  the  constant  hum  of  diesel  generators.  These  thirsty  machines  provide  the  base’s 
electric  lifeblood  enabling  lighting,  refrigeration,  life-support  heating  and  cooling,  battery 
charging,  and  equipment  operation.  Years  of  conflict  proved  our  dependence  on  fossil  fuels 
comes  with  a  cost.  A  constant  consumption  of  diesel  fuel  requires  a  hefty  logistics  tail 
to  provide  fossil  fuels  to  every  base  across  the  area  of  operations.  During  Operation  Iraqi 
Freedom  and  Operation  Enduring  Freedom  an  Army  study  found  that,  for  every  24  convoys 
sent  out,  one  service  member  or  civilian  engaged  in  fuel  transport  was  killed.  Additionally, 
due  to  transportation  cost  of  fossil  fuels  to  isolated  bases,  the  price  of  fuel  escalates  from 
approximately  one  dollar  per  gallon  to,  in  some  cases,  400  dollars  per  gallon  [1]. 

The  source  of  electricity  generation  service  members  depend  on  decreases  combat  capabil¬ 
ity.  Military  forces  are  taken  away  from  combat  operations  and  nation  building  to  guard 
and  operate  convoys.  Dollars  spent  on  expensive  fossil  fuels  are  better  allocated  to  support 
service  members  and  operations.  Additionally,  expeditionary  forces  do  not  have  the  ability 
to  transport  heavy  fuels.  Their  dependence  on  fossil  fuels  for  power  generation  inhibits 
their  ability  to  conduct  operations. 

This  paper  attempts  to  improve  upon  an  alternative  method  of  producing  electricity  while 
in  austere  conditions  without  the  necessity  of  fossils  fuels  and  their  accompanying  logistics 
tail. 
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Figure  1.1:  Soldier  guarding  fuel  convoy  after  attack  [2]. 


1.2  Piezoelectricity  Definition  and  History 


The  name  “piezo”  is  derived  from  the  Greek  piezo  or  piezein,  meaning  “to  squeeze  or 
press”  [3].  Piezoelectricity  is  pressure-driven  electricity.  Certain  materials  generate  elec¬ 
tricity  when  under  pressure  (mechanical  stress).  These  same  materials  change  shape  when 
electricity  is  applied  to  the  material.  A  common  use  of  the  piezoelectricity  phenomenon  is 
demonstrated  in  quartz  watches.  A  battery  in  the  watch  feeds  electricity  to  a  piezoelectric 
quartz  crystal.  The  electricity  induces  a  precise  vibration  in  the  quartz.  The  mechanisms  in 
the  watch  then  convert  the  quartz  vibrations  into  a  means  of  keeping  accurate  time. 
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Figure  1.2:  Piezoelectric  quartz  crystal  resonator  used  for  time  keeping  [4], 


In  1880,  two  French  physicists,  Pierre  and  Jacques  Curie,  found  that  in  certain  materials, 
such  as  zinc  blende,  sodium  chlorate,  boracite,  tourmaline,  cane  sugar,  Rochelle  salt,  and 
quartz,  mechanical  stresses  induce  macroscopic  polarization  and  hence  the  production  of 
electric  surface  charges.  The  following  year,  Gabriel  Lippmann  predicted  the  converse 
effect  from  fundamental  thermodynamic  principles,  that  is,  an  imposed  voltage  produces 
mechanical  deformations  or  strains  of  the  material  [5].  In  1882,  the  Curie  brothers  con¬ 
firmed  the  existence  of  the  converse  effect  based  on  experimental  observations  [3]. 

Pierre  Curie  used  the  piezoelectric  effect  to  measure  charges  emitted  by  radium,  but  piezo¬ 
electricity  was  not  used  for  practical  applications  for  several  decades.  During  World  War  I, 
Paul  Langevin  developed  piezoelectric  crystal  quartz  transducer  depth  sounding  devices  to 
locate  submerged  vessels,  primarily  German  submarines  [6].  In  the  1920s,  developments 
in  crystal  resonators  for  the  stabilization  of  oscillators  launched  the  field  of  frequency  of 
control.  The  introduction  of  quartz  control  drastically  changed  the  way  humans  keep  time. 
Today,  piezoelectric  applications  include  smart  materials  for  vibration  control,  aerospace 
and  astronautical  applications  of  flexible  surfaces  and  structures,  sensors  for  robotic  appli¬ 
cations,  and  novel  applications  for  vibration  reduction  in  sports  equipment  (tennis  racquets 


and  snowboards)  [5].  Recently,  scientists  have  focused  research  on  using  piezoelectric 
materials  for  energy  generation  [7]— [14] . 


1.3  How  Piezoelectricity  Works 

The  separation  of  positive  and  negative  electrical  charges  in  many  molecules  results  in 
what  is  known  as  a  dipole  moment.  In  piezoelectric  crystals,  when  a  mechanical  stress  is 
applied  (the  crystal  being  compressed,  twisted,  or  pulled)  the  molecular  dipole  moments 
re-orient  themselves.  This  re-orientation  causes  a  variation  in  surface  charge  density  and 
thus  produces  voltage.  The  effect  is  illustrated  in  Figure  1.3  for  forces  normal  to  the  mate¬ 
rial.  Conversely,  when  an  electric  field  is  applied  across  a  piezoelectric  medium,  there  is  a 
slight  change  in  the  shape  of  the  dipoles  causing  a  very  small,  but  significant  change  in  the 
material  dimensions. 

Pulling  normal  force  Equilibrium  Compressional  normal  force 

AAA 


V  V  V 

Figure  1.3:  Illustration  of  the  piezoelectric  concept  [3]. 


4 


1.4  Applications  of  Piezoelectric  Power  Generation 

Using  the  converse  piezoelectric  effect,  it  is  possible  to  generate  electricity  by  deforming 
piezoelectric  materials.  However,  piezoelectric  generation  is  far  from  the  most  efficient 
method  of  power  generation.  Presently,  piezoelectric  generators  can  only  produce  small 
amounts  of  electricity,  but  with  the  increase  of  efficiency  of  powered  devices,  piezoelectric 
power  generation  has  become  a  possible  supply  of  electricity  for  certain  devices. 

Devices  well  suited  for  piezoelectric  power  supply  include  the  following: 

•  devices  with  low  power  requirements 

•  devices  without  power  supplies  in  close  proximity 

•  devices  without  easy  access  to  a  power  supply  (power  lines  are  inconvenient  or  im¬ 
possible) 

•  devices  without  access  to  other  alternative  power  sources  (solar,  wind,  geothermal) 

•  devices  with  access  to  vibrations  or  other  mechanical  manipulation 

The  most  obvious  device  that  meets  the  above  criteria  is  a  sensor.  Sensors  are  often  located 
far  from  power  sources,  concealed  from  light  sources,  and  require  power  to  collect  and 
send  information.  It  is  not  always  possible  or  convenient  to  connect  power  to  the  sensors 
by  wire.  Examples  of  sensors  that  could  receive  power  through  piezoelectric  generation 
include  sensors  in  air  ducts  [15],  tires,  and  inside  machinery. 

Piezoelectric  power  generation  has  the  potential  to  enhance  military  operations  as  well. 
Military  units  often  operate  in  primitive/isolated  environments  where  it  is  not  practical  to 
transport  conventional  fuel  to  generate  electricity.  Military  units  must  rely  on  alternative 
methods  of  power  generation  or  do  without.  In  addition  to  solar  and  wind  power  genera¬ 
tion,  service  members  could  generate  power  with  piezoelectric  power  generators  fixed  to 
personnel,  equipment,  or  any  other  item  in  motion. 


1.5  Research  Objectives 

This  paper  attempts  to  find  the  optimal  designs  involving  material  parameters  in  the  piezo¬ 
electric  generator,  as  well  as  the  generator’s  environment,  in  order  to  maximize  electric 
output. 
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1.6  Thesis  Organization 

This  chapter  introduced  the  motivation  and  relevant  background  of  harvesting  piezoelectric 
energy.  In  Chapter  2,  we  review  a  one-dimensional  power  generation  model  and  study  how 
to  optimize  power  generation  using  variables  in  the  model.  From  the  results  of  Chapter  2 
we  derive  a  three-variable  power  generation  model.  In  Chapter  3,  we  derive  a  six-variable 
power  generation  model  using  a  single  mode  model.  In  Chapter  4,  we  describe  the  various 
methods  we  use  in  order  to  optimize  our  three-variable  and  six-variable  power  generation 
models.  In  Chapter  5,  we  discuss  the  results  of  our  optimization  methods  and  the  robust¬ 
ness  of  our  results.  In  Chapter  6,  we  conclude  the  paper  and  discuss  further  studies  to  be 
completed  based  on  the  information  and  ideas  generated  in  this  paper.  Finally,  we  include 
an  appendix  to  display  MATLAB  code  used  for  calculations  within  the  paper. 
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CHAPTER  2: 

One-Dimensional  Power  Generation  Model: 
Derivation  of  the  Three  Variable  Model 


2.1  1-D  Model  Design:  Interpretation  Scheme 

We  can  use  a  basic  1-D  model  to  analyze  the  power  generated  from  a  vibration  energy 
harvester  to  understand  conversion.  This  model  is  limited  to  harvesters  where  the  electri¬ 
cal  damping  term  is  linear  and  proportional  to  velocity.  Nevertheless,  this  simple  model 
is  useful  in  understanding  the  feasibility  of  the  device  and  input  parameters  on  power  ex¬ 
tracted.  The  electrical  energy  is  extracted  from  the  mechanical  system,  which  is  excited 
by  a  mechanical  input.  This  extraction  is  not  necessarily  linear,  or  proportional  to  velocity, 
however,  it  is  a  dissipative  process  and  can  generally  be  viewed  as  electrical  damping. 

Following  [16],  we  consider  a  generator  (shown  in  Figure  2.1)  which  consists  of  a  seismic 
mass  m  on  a  spring  k.  When  the  generator  is  vibrated,  the  mass  moves  out  of  phase  with  the 
generator  housing,  so  that  there  is  a  net  movement  between  the  mass  and  the  housing.  This 
relative  displacement  is  sinusoidal  in  amplitude,  and  can  drive  a  suitable  transducer  to  gen¬ 
erate  electrical  energy.  The  transducer  is  depicted  as  a  dashpot,  d ,  because  the  conversion 
of  mechanical  energy  into  electrical  energy  damps  the  mass.  There  are  several  transduction 
methods  suitable  for  the  generator.  The  choice  of  transducer  makes  little  difference  to  the 
amount  of  electrical  power  that  will  be  generated.  Three  possible  transduction  mechanisms 
are  (1)  piezoelectric:  using  piezoelectric  material  to  convert  strain  in  the  spring  into  elec¬ 
tricity;  (2)  electromagnetic:  a  magnet  attached  to  the  mass  induces  a  voltage  in  a  coil  as  it 
moves;  and  (3)  electrostatic:  an  electric  arrangement  with  a  permanent  charge  embedded 
in  the  mass  induces  a  voltage  on  the  plates  of  a  capacitor  as  it  moves. 
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Figure  2.1:  A  schematic  diagram  of  the  generator  [17]. 


In  the  simplified  model,  we  assume  that  the  mass  of  the  vibration  source  is  much  bigger 
than  the  mass  of  the  seismic  mass  in  the  generator,  and  the  vibration  source  is  an  infinite 
source  of  power.  This  implies  that  the  vibration  source  is  unaffected  by  the  movement  of 
the  generator. 

The  differential  equation  of  motion  can  be  derived  from  the  dynamic  forces  on  the  mass: 

mz  ( t )  +  dz  ( t )  +  kz  ( t )  =  -my"  ( t )  (2.1) 

where  the  generator  housing  is  vibrated  with  a  displacement  y(t),  the  relative  motion  of  the 
mass  with  respect  to  the  housing  is  z(t),  m  is  the  seismic  mass,  d  is  the  damping  constant, 
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and  k  is  the  spring  constant.  The  solution  to  this  second-order  ordinary  differential  equation 
(ODE)  with  constant  coefficients  consists  of  two  parts,  the  complementary  function,  which 
is  the  solution  of  the  corresponding  homogeneous  equation,  and  the  particular  solution.  The 
complementary  function,  in  this  case,  is  a  damped  free  vibration  which  goes  to  zero  as  time 
increases.  The  particular  solution  is  the  one  of  interest  here.  For  a  sinusoidal  excitational 
vibration,  y(t)  =  Yocos(cot),  the  particular  solution  to  the  preceding  ODE  is  a  steady-state 
oscillation  of  the  same  frequency,  ft),  as  that  of  the  excitation.  We  can  assume  the  particular 
solution  to  be  of  the  form 

z(t)  —  Acos(cot  —  <j>)  (2.2) 

where  A  is  the  amplitude  of  oscillation  and  (j)  is  the  phase  of  the  displacement  with  respect 
to  the  exciting  force. 

Differentiating  Equation  (2.2)  yields 


z'  (t)  =  —Aft)  sin  (tut  —  (/)) , 
z"  (t)  =  —Aft)2 cos  (cot  —  cj>) . 

Substituting  the  above  equations  into  Equation  (2.1)  we  have 


(2.3) 

(2.4) 


—  Amcor  cos  (cot  —  (j))  —Ad(0  sin  (cot  —  (j))  +A/:cos  (cot  —  (j>)=  mYoCO2  cos  (cot) .  (2.5) 


Utilizing  the  trigonometric  identities  cos  (cot  —  (j))  =  cos(C0t)cos((j>)  +  sin((Dt)sin(ty)  and 
sin(oot  —  <f>)  =  sin(cot)cos(<f>)  —  cos((Ot)sin(<f>)  we  can  rewrite  Equation  (2.5)  as 


— Am  co2  cos  cot  cos  (f)  —  AmCO2  sin  cot  sin  0  —  AJ ft)  sin  cot  cos  0  +  Ad  cocos  (Ot  sin  (f) 
+Ak  cos  cot  cos  (j)  +Aksincot  sin^  =  mY^co2  cos  (cot) . 


(2.6) 


Matching  the  coefficients  for  cos  (cot)  and  sin  (cot),  respectively,  we  find 


(Ak  —  Amco2)  cos  (j)  +A<ift)sin0  =  mYoCO2, 

(2.7) 

(k  —  mco2)  sin  (j)  —  dcocoscj)  =  0. 

(2.8) 
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This  leads  to 


tan  (j)  = 


dco 


k  —  m(0 


2’ 


[  (k  —  mco2)  cos  (j)  +  dco  sin  (j))  A  —  mYoCO2 , 


[  (k  —  mco2)  +  dco  tan  cf>]  A  =  mYoCO 2  sec  (j)  =  mYoCO 2  \J  1  +  tan2^ . 
Solving  for  A,  we  obtain 


(2.9) 

(2.10) 

(2.11) 


inYo  co2  a/1  +  tan20 
(k  —  mco2)  +dcotan(f) 


mYoCO2 _ 
{k-m(02)+d(o(j^j 

mY()(02 

A  —  ,  = 

y  (k  —  mco2)2  +  (dco)2 


mYoCO2  \J  (k  —  mco2)2  +  (dcoY 
(k  —  mco2)2  +  (dco)2 


(2.12) 

(2.13) 


We  now  express  the  amplitude  and  phase  in  non-dimensional  form.  Dividing  the  numerator 
and  denominator  by  k,  we  get 


tan0 


dco 

k  —  mco2 


dco 

k 


m(Q 2  5 
k 


(2.14) 


A  = 


(2.15) 


These  expressions  can  be  further  expressed  in  terms  of  the  following  quantities  where  COn 
is  the  natural  frequency  of  undamped  oscillation  (resonant  angular  frequency),  dc  is  critical 
damping,  and  £  is  the  modal  damping  ratio: 


CO 1 1  — 


(2.16) 


dc  —  IniCO,,. 


dco 

k 


d  dcCO 
dc  k 


2C’ 

COn 


(2.17) 

(2.18) 
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(2.19) 


k  k 
dc  2m  ft),. 


The  non-dimensional  formulas  for  the  amplitude  and  phase  then  become 


k/m  I  k  ft),, 

0  [~k  2  V  m  2 

i\  — 


2t  — 

tan  cj)  = - 


1-  -® 
V  co„ 


(2.20) 


(2.21) 


The  instantaneous  power  transfer  in  the  mass  is  the  product  of  the  force  on  the  mass  and  its 
velocity: 

P  (t)  =  -my"  (; t )  [ y '  (. t )  +  z  (f )]  .  (2.22) 

When  damping  is  present,  due  to  the  electrical  transducer,  there  is  a  net  transfer  of  mechan¬ 
ical  power  into  electrical  power.  This  net  electrical  power  generated,  P,  is 

T  T 

P  =  —  f  p(t)dt  —  —  j  mco2Yo cos  (cot)  [—YqCO  sin  (cot) —Aco sin  (cot  —  cj))]dt,  (2.23) 


1  4-  cos  (2  cot) 


P—  —  J  mAco  Yo sin (f) cos  (cot)dt  —  —  j  mAco  Yo  sin  ({> - — 

0  0 


dt ,  (2.24) 


P  =  mYoCO3  —  sin  cj). 

Substituting  in  Equation  (2.21)  for  A,  we  obtain: 


(2.25) 


P  =  mY02CO3- 


2V  ‘-5  + 


2  y/l  +COt 2(j) 


(2.26) 
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P  =  mY02(o3- 


1-5)  +(2^ 


(2.27) 


l-f 

\(0n  ) 

2C  — 


P  = 


mYo2CD 3 


CO  \  CO 

(On  I  fix. 


ft>- 


1 


5)  +  (2^ 


1 


5)  +  (2^ 


2  ‘ 


(2.28) 


We  now  have  an  explicit  expression  for  net  electric  power  generated  we  can  conduct  analy¬ 
sis  on.  In  the  following  sections  we  first  view  the  output  power  as  a  function  of  normalized 
frequency  and  investigate  how  the  output  power  behaves  as  we  vary  normalized  frequency 
for  various  values  of  the  modal  damping  ratio.  After  that,  we  treat  the  output  power  as 
a  multi-variable  function  of  both  normalized  frequency  and  modal  damping  ratio,  seeking 
the  optimal  value  of  the  output  power.  Then,  we  compare  our  model  prediction  with  an 
experimental  result. 


2.2  Power  as  a  Function  of  Normalized  Frequency 

First,  we  will  treat  net  electric  power  generated,  P,  as  a  function  of  normalized  frequency 
only.  To  do  this,  we  will  introduce  C  =  m^Y(jd)3  and  x  =  ^  to  get  the  function 


P(x) 


C 


(1  —x2)2  +  (2£x) 


T 


(2.29) 


In  order  to  optimize  the  function  we  differentiate  P(x)  with  respect  to  v: 


2Cx5  [x4  —  4x2  +  8£2x2  +  3] 
(. x 4  -I-  4£2.v2  —  2.v2  +  l)2 


(2.30) 


By  inspection,  we  can  see  if  £  =  0,  P'(x)  does  not  exist  when  .r  =  1.  In  fact,  P(l)  becomes 
infinitely  large.  If  £  ^  0,  the  critical  points  of  P(x)  depend  on  the  value  of  £. 
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2.2.1  Initial  Analysis 

We  will  fix  values  for  £  and  find  the  value  of  x  that  maximizes  P(x).  Like  Williams  and 
Yates  [17],  we  will  analyze  P(x)  with  £  values  of  0.1, 0.3, 0.5,  and  0.7. 

For  £  =  0.1, 


2 Cx5  [x4-3.92x2±3] 
(x4  —  1.96x2  ±  1)" 


(2.31) 


Solving  P\x)  =  0  yields  critical  points  x  —  ±1.6963,  ±1.021 1.  Recall  that  x  =  ^  >  0.  For 
£  =  0.1,  we  plot  (^(2  as  a  function  of  x  in  an  interval  containing  the  critical  points  in  Figure 
2.2. 


From  Figure  2.2,  we  can  see  that  the  maximum  power  can  be  generated  when  x  =  = 

1.0211  for  £  =  0.1,  and  the  maximum  value  of  is  26.0421. 
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For  C  =  0.3, 


.  2Cx5\x4- 3.28.x2 +  3] 

P  W  =  - k - T1- 

(. x 4  —  1.64x2  +  1) 

Solving  P'(x)  —  0  yields  four  complex  roots.  The  critical  point  is  x 
for  C,  —  0.3  in  Figure  2.3. 


(2.32) 

0.  Now  we  plot 


Figure  2.3  shows  when  £  =  0.3,  is  an  increasing  function  of  x. 
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For  C  =  0.5, 


P\x) 


2Cx5  [ x 4  —  2x2  +  3] 
(x4  —  x2  + 1)“ 


(2.33) 


Solving  P'(x)  —  0  yields  four  complex  roots.  The  critical  point  is  x  =  0.  Now  we  plot 
for  £  =  0.5  in  Figure  2.4. 


Figure  2.4  shows  when  £  =  0.5,  is  an  increasing  function  of  x. 
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For  C  =  0.7, 


.  2C.r5  \x4  —  0.08x2  +  3] 

p'  (  v)  —  _ h _ I 

(*4-0.04x2  +  l)2 

Solving  P'(x)  —  0  yields  four  complex  roots.  The  critical  point  is  x 
for  £  =  0.7  in  Figure  2.5 


(2.34) 

0.  Now  we  plot 


Figure  2.5  shows  when  £  =  0.7,  js  an  increasing  function  of  x. 
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2.2.2  Analysis  as  £  Approaches  Zero 

One  will  notice  only  £  =  0.1  resulted  in  a  critical  point  other  than  x  —  0.  Let  us  evaluate 
the  critical  points  and  corresponding  optical  when  £  <0.1.  For  £  =  0.05, 


2Cx5  [x4-3.9Sx2  +  3] 
(x4-1.99x2  +  lf 


(2.35) 


Solving  P'(x)  —  0  yields  critical  points  x  —  ±1.7233,3=1.0051.  For  £  =  0.05,  we  plot 
as  a  function  of  x  in  an  interval  containing  the  critical  points  in  Figure  2.6. 


From  Figure  2.6,  we  can  see  that  the  maximum  power  can  be  generated  when  x  =  = 

1.0051  for  £  =  0.05,  and  the  maximum  value  of  is  101.01. 
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It  appears  as  £  approaches  zero,  the  optimal  value  of  x  approaches  one  and  the  resulting 

P(x) 

maximum  power  generation,  increases.  See  Figures  2.7  and  2.8  for  further  illustration 
of  the  how  the  optimal  x  value  and  maximum  power  generation  change  as  £  decreases. 


Figure  2.7:  Optimal  x  as  a  function  of  £. 


Figure  2.7  further  depicts  the  optimal  frequency  ratio  (^-)  approaching  one  as  £  approaches 
zero.  This  follows  the  principal  of  resonance,  where  a  vibration  is  able  to  drive  a  sys¬ 
tem  into  larger  oscillations  when  the  vibration  frequency  matches  the  system’s  natural  fre¬ 
quency  (co  =  co,  J. 
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2500 


P(x) 

Figure  2.8  demonstrates  that  when  we  set  a  =  C0,„  -j1  depends  solely  on  how  small  we 
can  set  £. 

2.3  Power  as  a  Function  of  Modal  Damping  Ratio  and 
Normalized  Frequency 

We  will  now  treat  net  electric  power  generated,  P,  as  a  function  of  both  modal  damping 
ratio,  C,  and  frequency  ratio,  x  =  jy.  We  can  graphically  solve  for  the  optimal  value  of 
by  plotting  as  the  dependent  variable  and  using  £  and  x  —  as  the  independent 
variables.  If  we  use  the  parameters  0.1  <c<  1  and  0  <  x  <  2,  we  get  the  plot  in  Figure 
2.9. 
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From  Figure  2.9  we  can  see  the  extreme  value  of  increases  as  £  goes  to  zero  and  x  goes 
to  one. 


„  10  , 

C  X=0/ro 

^  n 


Figure  2.9:  Power  as  a  function  of  damping  factor  and  frequency  (0.1  <c<  1  and  0  <  *  <  2) 
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The  extreme  value  of  net  electric  power  generated  is  further  illustrated  in  Figure  2.10  with 
parameters  0.05  <  £  <  1  and  0  <  x  <  2. 


Figure  2.10:  Power  as  a  function  of  damping  factor  and  frequency  (0.05  <c<  1  and  0  <  x  <  2). 


2.4  Model  Comparison  with  Experimental  Results 

Roundy  and  Wright  [15]  use  similar  mathematical  methods  to  design  and  test  a  piezoelec¬ 
trical  vibration  based  generator.  Their  design  uses  an  8.44  gram  mass  and  a  two-layer  sheet 
of  PZT-5A  with  a  steel  center  shim  excited  by  vibrations  of  2.5 m  ,v  2  at  120 Hz-  They 
assume  the  resonance  frequency  of  their  generator,  con,  matches  the  driving  frequency,  co. 
The  damping  ratio,  £,  is  measured  by  applying  an  impulse  to  the  system,  and  then  measur¬ 
ing  the  open  circuit  voltage  output.  The  resulting  damped  harmonic  oscillation  is  used  to 
calculate  a  damping  ratio  near  0.015  [15].  Figure  2.11  shows  the  piezoelectric  generator 


21 


designed  by  Roundy  and  Wright. 


Figure  2.11:  Piezoelectric  generator  prototype  used  by  Roundy  and  Wright  [15]. 


With  these  parameter  values,  we  can  compute  the  power  generated  using  our  ODE  model. 
All  red  text  below  are  values  taken  from  the  Roundy  and  Wright  generator  model. 


(0  —  (on  —  {2k)\2QHz 


second 


Amplitude 


a 

x  —  —  =  1 

(On 


£  d  d 
dc  2  moon 


=  0.015 


m  «  8.44 grams  =  8.44  x  10  5  kg 
2.5  m 

acceleration  = 


>o  = 


acceleration 


(Oz 


second 2 

2.5  m 
second 2 

(  240 7l  y 
\  second  ) 


=  3.518  x  10' 


in 


(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 
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C  -  m^Ytcol  =  (8.44  x  lCT3kg)(0.015)(3.518  x  10“6m)2  f  (2.42) 

\second  J 


C  = 


1.0494  x  10-6%  ■  tn2 


(2.43) 


second 3 

Substituting  values  given  from  Roundy  and  Wright  and  values  calculated  above  into  Equa¬ 
tion  (2.29),  we  get: 


n  i)  = 


1.0494  x  10 ~6kg-m2 
second 3 


1 


1.2  x  \0~3kg-m2 


(2(.015))' 


second 3 


(2.44) 


The  corresponding  result  from  Roundy  and  Wright  is 

3.75  x  10 ~4kg-m2 


P(  1) 


second 3 


(2.45) 


Roundy  and  Wright’s  experiment  results  in  power  68.75%  less  than  our  mathematical  re¬ 
sults.  The  disparity  between  our  ODE  model  calculation  and  Roundy  and  Wright’s  exper¬ 
imental  results  can  be  explained  through  the  fact  our  model  is  purely  mechanically  based 
and  does  not  take  into  account  loss  of  energy  through  the  transfer  of  mechanical  energy 
to  electrical  energy.  Despite  the  power  disparity,  Roundy  and  Wright’s  results  support  our 
ODE  model. 


2.5  One-Dimensional  Power  Generation  Model  Conclu¬ 
sion 

Our  1-D  model  shows  potential  to  serve  as  a  foundation  to  optimize  power  generated  by 

p(x\ 

a  piezoelectric  generator.  In  order  to  optimize  power  generated,  -^4,  we  should  minimize 
the  damping  factor,  C,  and  choose  the  optimizing  frequency  ratio,  x  —  77-.  Keep  in  mind 
that  as  £  approaches  zero  the  optimal  frequency  ratio  approaches  one  and  the  maximum 
power  generated  drastically  increases. 
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CHAPTER  3: 
Single-Mode  Model: 
Derivation  of  the  Six  Variable  Model 


There  are  many  studies  on  the  transduction  of  electrical  energy  from  a  single  piezoelectric 
transducer  (see  [18]  and  references  therein)  and  the  general  conclusion  is  that  piezoelectric 
energy  harvesting  approaches  have  considerable  promise  for  applications  involving  small 
amounts  of  power  [19].  To  boost  the  power  output  of  such  system,  one  proposal  by  Scruggs 
[20]  is  to  use  multiple  transducer  patches  and  channel  energy  generated  by  some  of  the 
patches  into  other  patches  in  order  to  excite  additional  vibration  modes.  An  example  of 
this  type  of  energy  harvesting  system  is  depicted  in  Figure  3.1.  The  parameters  describing 
the  size  of  a  piezoelectric  patch  take  the  typical  values  l\  ~  I2  ~  h  —  100 mm,  w  —  25 mm, 
tp  =  0.25 mm  and  the  thickness  of  the  beam  is  t b  =  0.25 mm  [21]. 


t\ - 4* - ^ - 4* - - 4 


Figure  3.1:  Energy  harvesting  bimorph  cantilever  with  distributed  piezoelectric  transducers  (gray) 
bonded  to  a  substrate  (white).  Adapted  from  Dutoit,  Wardle  and  Kim  [16]. 
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The  general  equations  for  linear  piezoelectricity  are  shown  in  Figure  3.2. 


Gauss'  law  : 


cD 


~  =  P- 


CXi 

\D' 

D\x\=  D-, 

LA. 

p.  charge  densit>r 

Mechanical  energy  balance:  /; 

u. :  displacement.  T..  :  stress 

Constitutive  relation  [2]: 


:  a  vector  of  electrical  displacements  (charge  area) 


c  u. 


ct‘ 


ar. 


cx 


[A" 

^li 

t? 

O 

_ 1 

’A" 

:  A 

^11 

^5 

e2 

A 

e33 

*1 

«3l 

e33 

E3 

71 

~e3l 

C11 

C12 

C13 

sl 

r2 

= 

~e3l 

^1 

C13 

s2 

t3 

~e33 

C13 

c13 

c33 

s3 

7 ; 

~el5 

C55 

s4 

T \ 

“e15 

C55 

S5 

1 

K? 
_ 1 

.  0 

Ss 

s(x)=\ 


S1 

A~ 

'7]' 

s2 

t22 

t2 

= 

s3 

s4 

material  strains,  T 1  a*  1  = 

T33 

t23 

t3 

z \ 

S5 

A 

7] 

s6 

.A. 

material  stress  (force  area). 


E(x)  = 


electrical  field  in  the  material  (volts  meter) 


Figure  3.2:  General  equations  for  linear  piezoelectricity  [20]. 
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The  constitutive  relation  in  Figure  3.2  is  for  material  that  is  electrically  and  mechanically 
isotropic  in  the  1  and  2  directions,  while  the  3-direction  is  the  poling  direction.  The  first 
three  rows  of  the  constitutive  relation  relate  the  electrical  displacement  to  the  electrical  field 
and  to  the  strain  in  the  material,  and  the  remaining  six  rows  relate  the  mechanical  stress  to 
the  electric  field  and  to  the  strain.  Here  the  stress,  T(x),  and  strain,  S(jc),  are  represented  in 
Voigt  notation  or  Voigt  form.  Referring  to  the  Figure  3.1  above,  1  corresponds  to  x  and  3 
corresponds  to  y. 

Equations  for  a  patch  on  a  beam  are  largely  based  on  Euler-Bernoulli  (or  Timoshenko) 
beam  models.  Piezoelectric  layers  or  patches  are  bonded  to  the  both  sides  of  the  beam.  The 
beam  is  assumed  to  be  in  pure  bending;  all  other  deformations  are  considered  negligible. 
Here  we  consider  a  “bimorph”  configuration.  That  is,  in  the  undisturbed  state  the  beam  lies 
along  the  x-axis  in  a  <  x  <  b,—  ^  <  y  <  j  <  z  <  j,  with  equal  patches  on  its  upper 
and  lower  surfaces  y  =  ±^,  each  of  thickness  tp  and  covering  the  entire  upper  and  lower 
surfaces.  It  is  further  assumed  that  the  patches  are  bonded  perfectly  to  the  substructure 
and  the  electrodes  cover  their  entire  upper  and  lower  surface  areas.  The  beam  is  clamped, 
at  x  =  a,  on  a  base  that  is  subject  to  motion  in  the  y-direction;  it  is  free  at  x  =  b.  The 
polar  direction  on  the  top  patch  is  in  the  y-direction  and  on  the  bottom  patch  it  is  reversed 
so  that  the  voltages  on  the  two  patches  add.  Due  to  this  asymmetry,  we  can  compute  the 
piezoelectric  effect  on  the  upper  patch  alone  and  double  the  result. 

Adopting  the  standard  assumptions  of  Euler-Bernoulli  beam  theory,  a  balance  of  forces  and 
moments  can  be  combined  to  yield 

M-xx  —  m  t  T  htt)  ■  (3.1) 

where  M(x,  t )  is  the  internal  moment  generated  by  mechanical  and  electrical  strain,  m  is  the 
mass  per  unit  length  of  the  composite  beam,  co  is  the  relative  vertical  deflection  along  the 
beam  from  the  base,  and  h  is  the  absolute  motion  of  the  base/host  structure  (thus  w  +  h  is 
the  absolute  displacement  of  the  beam).  The  Euler-Bernoulli  bending  deformation  model 
relates  the  stress  to  the  curvature  of  the  beam: 


Si  =  -y(Oxx 


(3.2) 
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The  constitutive  relations  yield  7)  =  cnSi  Te3iE3,  D3  =  ±^31^1  +  e33£3  (the  sign  of 
the  coupling  terms  may  vary  from  patch  to  patch  depending  on  whether  they  are  poled 
in  the  positive  y  or  negative  y  direction).  Additionally,  it  is  assumed  that  the  electrical 
field  is  constant  throughout  the  thickness  of  the  patches;  therefore,  £3  =  ±p  with  the  sign 
depending  on  the  poling  direction,  where  vp  is  the  terminal  voltage  on  the  top  of  the  top 
patch  with  thickness  tp. 

The  relation  between  moment  and  stress  is: 

[b 
2 

M  =  —  w  I  T\ydy—w 

_'± 

2 

M  =  c  1 1  w  ~~  (Oxx  -  g3iw— + tp  [H(x-a)-H(x- b)]  vp  ■  2,  (3.4) 

where  w  is  the  width  of  both  the  beam  and  the  piezo  patches  and  H  denotes  the  Heaviside 
step  function. 

When  the  electric  field  is  normal  to  the  beam  axis  and  uniform  in  the  patch  over  a  <  x  < 
b  and  zero  outside  this  interval,  the  piezoelectric  effect  contribution  to  the  displacement 
equation  enters  as  derivatives  of  delta  functions  at  the  ends  of  each  patch  [22] : 

m(Ott+  (Oxxxx  +  e3]w(th  +  tp)  [d'(x-b)  -  8'(x-a)\  vp  =  -mhtt.  (3.5) 


The  boundary  conditions  for  equation  (3.5)  are  ( o(a,t )  =  CQx(a,t )  =  0.  This  reflects  the  fact 
that  the  beam  is  clamped  at  x  —  a,  but  at  the  free  end,  x  =  b,  we  have  M  =  Mx  —  0.  The 
size  of  the  charge  can  be  found  by  integrating  the  Gauss’  law  over  the  surface  area  of  the 
electrodes  [22]: 

q  (t)  =  If  D3dA-  ff  D3dA  =  e3\w(tp  +  ts)  [d)x\  -  — (b  -  a)  vp  .  (3.6) 

JJ  upper  JJ  lower  a  tp 

Equations  (3.5)  and  (3.6)  provide  a  complete  system  of  equations  for  the  vibration  energy 
harvester.  The  most  common  method  of  solving  the  system  (3. 5-3. 6)  is  to  assume  that  the 
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beam  deflection  can  be  represented  by  an  eigenfunction  expansion 

oo 

ft)  (x.t)  =  Y,  0/  W  Mi  (0)  (3-7) 

/=1 

where  fa{x)  is  the  i-th  mode  shape  function  and  «,•(/)  is  the  i-th  modal  displacement.  The 
mode  shape  functions  must  satisfy  the  boundary  conditions  and  follow  the  following  or¬ 
thogonality  relationships: 

b 

I  (f)i  (x)  <j>j  (x)  dx  =  8jj.  (3.8) 

a 

Substituting  (3.7)  into  (3.5),  multiplying  both  sides  by  fa(x),  and  integrating  x  from  a  to  b 
yields  for  each  k,  we  arrive  at 


Uk"  0)  +  [wscm]  uk  +Akvp  =  -m  [■ yk }  h"  (t) ,  (3.9) 


where  Equation  (3.10)  is  the  modal  short-circuit  frequency,  Equation  (3.11)  is  the  modal 
electromechanical  coupling  coefficient,  and  Equation  (3.12)  is  the  modal  influence  coeffi¬ 
cient  of  the  base  excitation. 


wSC,k 


V 


Cii  wtb3 
12  m(b  —  a) 


(3.10) 


where  Xk  is  the  eigenvalue  corresponding  to  (j)k(x);  for  cantilevered  configurations,  these 
numbers  are  the  solutions  to  cos(A)  cosh(A)  =  1 


Ak  =  c3iw  {tb  +  tp )  [$k  {x-b)~  fa'  {x  -  a)] ,  (3.11) 

b 

Yk  =  J  fa(x)dx.  (3.12) 

a 

At  this  point,  it  is  common  to  add  in  a  modal  damping  term  2 wsc,kfauk  (?)  to  the  left  side 
of  Equation  (3.9)  to  account  for  all  proportional  damping  effects.  Four  possible  damping 
mechanisms,  one  external  and  three  internal,  and  various  combinations  of  these  mecha¬ 
nisms  have  been  considered  by  Banks  and  Inman  [23].  They  are  viscous  air  damping,  strain 
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rate  damping  or  Kelvin- Voigt  damping,  spatial  hysteresis,  and  time  hysteresis,  respectively. 
Their  study  shows  that  air  damping  plays  a  more  significant  role  in  lower  modes  whereas 
internal  damping  plays  a  more  important  role  for  higher  modes. 

Equations  (3.9)  and  (3.6)  provide  a  complete  system  of  equations  for  the  vibration  en¬ 
ergy  harvester.  If  the  beam  is  excited  with  a  driving  frequency  that  is  close  to  one  of 
its  natural  frequencies,  then  the  corresponding  mode  of  that  frequency  will  dominate 
the  motion  of  the  beam.  In  this  situation,  it  is  reasonable  to  make  the  approximation 

oo 

(o(x,t)  —  £  (f)j(x)ui(t)  ~  0i  ( x)u\  ( ? ).  After  we  drop  the  ”1”  subscript,  this  simple  ap- 

i=  I 

proximation  reduces  Equation  3.9  to 

u"  (?)  +  2wSc^u  (?)  -1-  [wSc2]  u  +Avp  =  -m  [y]  h"  (?) .  (3.13) 


Here  m  is  the  mass/unit  length,  y  is  the  modal  influence  coefficient  of  the  base  excitation, 
£  is  the  modal  damping  ratio,  wsc  is  the  modal  short-circuit  natural  frequency,  and  A  is  the 
electromechanical  coupling  coefficient  (effective  piezoelectric  constant).  Equation  (3.6)  is 
simplified  to 


<?(?)  =£>31  W(tp  +  ts)  [<M 


2C33 


w(b  —  a)  v 


p  • 


(3.14) 


Differentiating  Equation  (3.14)  yields 


—  Au  4 - —  w  ( b  —  a)  Vp  —  —i. 


(3.15) 


where  i  represents  the  current.  We  can  focus  on  the  single-mode  model  (Equations  (3.13) 
and  (3.15))  for  now,  which  assumes  that  the  beam  is  vibrating  near  its  fundamental  natural 
frequency  and,  consequently,  the  motion  of  the  beam  is  described  by  one  modal  coordinate 
(typically,  the  fundamental  mode).  We  need  to  solve  the  system  and  find  the  maximum 
power  points  of  energy  harvesting  over  a  range  of  base  excitation  frequencies,  patch  thick¬ 
ness  and  length,  load  resistances,  etc.  The  results  will  allow  a  designer  to  choose  the 
optimal  resonant  frequency  and  patch  dimensions  to  maximize  the  power  harvested. 


In  the  single-mode  model,  i.e.,  Equations  (3.13)  and  (3.15),  the  two  unknowns  are  the 
modal  displacement  u(t)  and  the  terminal  voltage  vp,  and  the  two  inputs  are  the  base  motion 
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h{t)  and  the  current  i.  The  beam  tip  deflection,  which  is  a  predominant  measure,  can 
be  computed  by  the  simple  relation  c o{L,t )  =  (j)  { L)u(t ),  where  (p{x)  is  the  fundamental 
mode  shape  function.  The  coefficients  of  Equations  (3.13)  and  (3.15)  can  be  calculated 
using  knowledge  of  the  materials  and  the  layout  of  the  beam,  or  they  may  be  estimated 
experimentally. 


The  most  commonly  analyzed  scenario  for  vibration  energy  harvesting  is  steady  state  base 
excitation.  At  steady  state,  we  suppose  that  the  beam  and  base  motions  are  of  the  following 
forms,  respectively 

u(t)=U  sin  {(Obaset) ,  (3.16) 

h  it)  H  sin  {(Obaset  (frbase )  j  (3.1V) 

where  U  is  the  amplitude  of  the  modal  displacement,  H  is  the  amplitude  of  the  base  motion, 
(Obase  is  the  base  motion  driving  frequency,  and  (phase  is  the  phase  lead  of  the  base  motion 
relative  to  the  beam  motion.  If  we  further  consider  a  simple  resistive  load  where  i  —  ^  and 
R  is  the  load  resistance,  then  Equations  (3.13)  and  (3.15)  are  simplified  to  a  linear  system. 
Assuming  the  voltage  signal  takes  the  form  vp  =  Vp  sin ((%,„./  —  6voitage)  and  substituting 
these  expressions  into  Equation  (3.15),  we  have 


AU  (Obase  COS  ( (Obaset )  3“  CqV p(Obase  COS  ( (Oj}  aset  ® voltage ) 


Vp  sin  {(Obaset  @ voltage ) 

R 


(3.18) 

where  Co  =  f33vr  {b  — a).  With  the  help  of  trigonometric  identities  cos  (x  —  y)  =  cos* cosy  + 

tp 

sin*siny  and  sin(*  — y)  =  sinvcosy  —  cos*siny  we  can  rewrite  Equation  (3.18)  as 


AU  (Obase  COS  {(Obaset)  A  C()Vp(Obase  [cOS  {(Obaset)  CC)S  dygltage  T  sin  {(Obaset)  sin0vo/fflgeJ 
=  -j^  [sin  {(Obaset)  COS  Qyoltage  COS  {(Obaset)  sin  QVoltage\  • 

(3.19) 

Matching  the  coefficients  of  cos  {(Obaset)  and  sin  {(Obaset)  respectively,  we  obtain 

Vp 

AU  (Obase  +  CqVp  (Obase  COS  dvoltage  =  sin  Ovo/uige-  (3.20) 

Solving  for  cos  0voitage,  we  obtain 


COS  0voltage  —  C()(OixlseR  sin  QVoluige- 


(3.21) 
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Substituting  Equation  (3.21)  into  Equation  (3.20),  we  get 


-  AU (Bbase  3“  C()Vp(l)base  (  C()0)fjaseR  0  voltage)  ^  sin  QVoltage-  (3.22) 


Solving  for  sin  9voitage,  we  obtain 


' voltage 


AU  (ObaseR 

Vp(l+C02cobase2R2)' 


Substituting  Equation  (3.23)  into  Equation  (3.21)  yields 


(3.23) 


' voltage  — 


CpAU  CObase^R2 

Vp  (l  +  Cq2 (Obase2R2) 


Now  substituting  Equations  (3.23)  and  (3.24)  into  Equation  (3.13) 

U  (Obase  sin  (  (ObaseR  3”  2w SC  C  U  (Obase  COS  (  (ObaseR 
3”  [wsc]  U  sin  {(Qbaset)  3-AVpSin  (( Obaset  @ voltage ) 
=  /ft  [y]  H (Obase  sin  ( (Obaset  tybase)  i 


(3.24) 


(3.25) 


which  leads  to 

U (Obase  sin  (( Obaset)  3“  2wg(gl^U (Obase  COS  {(ObaseR  3“  [wsC-]  U  sin  (( Obaset ) 
3 -AVp  [sin  ((Obaset)  COS  Q voitage  -  COS  (ft)/,  aset )  sill  9voltage\ 

/ft  [y]  H (Obase"  [sift  ( i Obaset )  COS  < Phase  CO S  ( ( Obaset')  sift  tybase]  ■ 

Matching  the  coefficients  of  cos  ((Obaset)  and  sin  ((Obaset),  respectively,  we  get 


(3.26) 


'U  (Obase  3“  [ft\SC  ]  3~  AVp  COS  0  volt  age  /ft  [y]  H  (Obase  ^Q^tybase: 


(3.27) 


2wscCU(Obase  AV p  sill  0  volt  age —  bbl  [y]  H  (Obase  si  n  (f>ba 

Substituting  Equation  (3.21)  into  Equation  (3.27),  we  obtain 


(3.28) 


U  (Oba 


CqAU  (Obase"  R" 


:  .  r  21  TT  I  11/  '-'U vjbase  “ 

+  [wsc  u  +AVPVp  (1  +C„2ffl,„„2R2) 


=  m  [y\H  (Obase2  C  os  <j>ba 


(3.29) 
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Solving  for  cos  (j)base,  we  find 


2  |  1"-.*,  2"|  |  CqA  C0[yase  R 

-«w+KI  +  (iss;sij 

COS  <j)base  = - —— - « - U. 

Ill  [y]  H  CQbase" 

Substituting  Equation  (3.23)  into  Equation  (3.28),  we  obtain 

2w’sc  c  U  (Obase  +  AVp  AU®>™eR  =  _m  [y]  H  CQbase2  sin  (j)ba 

Vp  [l+Co  cobase  R  ) 

Solving  for  sin  (j)base,  we  find 


(3.30) 


(3.31) 


2wsc^(oh 

ase  + 


A  ®basc’R 


sin  <j)ba 


(l+Q)2  (Obase2  R2) 


m  [y]  H  (Oba 


U. 


(3.32) 


Using  Equations  (3.30)  and  (3.32)  and  the  trigonometric  identity  cos2 ((j)base)  +  sin2 ((j)base)  — 
1,  we  find 


m,  2  4-  L,  21  4-  CQA2Wbasl,2R 2 
CQbase  +  K-C  j  +  (t+0,2^2^) 

/n  [y]  Hcobase~ 


1  2 


u2  + 


+ 


^  ®baseR 


(  l+Co^Wiase2^2) 


m  [y]  H  CQha 


U2=l. 


(3.33) 


Solving  for  U  leads  to 


jj  w[y\CC 0)j)ase  (l +C(j  COj)ase  7?“) 

\/(  (  ^Ibasf-  '  [w£C-]  )  (  1  +Co"  CObase~R~  )  l  C()A  -  (Obase^R~  )  I  (-M'ScC  (Obase  (  1  (Obase~  R~  )  -\-A~  (QbaseR^  ~ 

(3.34) 

Using  Equations  (3.23)  and  (3.24)  and  the  trigonometric  identity  cos2(6voitage)  +  sm2(6voitage) 
1,  we  obtain 


AU  (QbaseR 

2 

1 

CqAU  CQbase2  R2 

Vp  (l  +  Co2  (Obase2 R2) 

i 

Vp  (l  +  Co2  CQbase2 R2) 

(3.35) 


Solving  for  Vp,  we  find 


Vp  —  AU  (Qbasc,R. 


(3.36) 
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Substituting  Equation  (3.34)  into  Equation  (3.36),  we  obtain 


Rtti  y  H  0)j,us,  (l+Co  (Ofrase~R 


((- ®base2+  WSC2])  {i+Co2 CObase2 R2)+CqA2 CObase2 R2)~ +  (2wsc£®base(l+Co'> CQbase2R2)+A2CObaseRy 


(3.37) 


The  average  power  dissipated  by  the  load  resistor  is  given  by 


V2 

P=-P~. 
2  R 


(3.38) 


Substituting  Equation  (3.37)  into  Equation  (3.38),  we  obtain  an  expression  for  power: 

p_  _ A2Rm2 [y]2H2 cobase6 ( 1 +C02 cobase2R2 ) 2 _ 

“  [(  (  ^bas(2  ^  [^;:2]  )  (  1  rC()~  ~  C{\A- 0)j}/!slr  ^  (2u'.vcC  ®base(  1  I  Q )  ~  G^base  ""  R 2  )  "t  ^  “  ^bas  e  ^)  ~  J 

(3.39) 

Equation  (3.39)  provides  an  explicit  and  complicated  form  for  power  generation,  which 
involves  nine  parameters.  We  will  study  the  effects  of  these  parameters  in  the  following 
chapters. 
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CHAPTER  4: 

Optimization  Methodology 


4.1  Optimization  Methodology  Overview 

We  optimize  piezoelectric  power  generation  by  analyzing  two  power  generation  models 
using  two  sampling  approaches  and  one  iterative  approach.  The  goal  of  our  methodology 
is  to  discover  the  optimal  values  of  each  variable  in  order  to  maximize  power  generation. 
In  the  following  chapter  we  compare  the  results  of  each  optimization  method  and  analyze 
the  robustness  of  our  results. 


4.2  Power  Generation  Models 

We  optimize  piezoelectric  power  generation  by  optimizing  the  mathematical  models  de¬ 
scribed  in  Chapters  2  and  3.  We  attempt  to  find  the  optimal  combination  of  variable  values 
in  order  to  maximize  power. 


The  first  model  is  an  adaptation  of  Equation  (2.28).  It  consists  of  three  variables,  ft),  ft),, , 
and  £.  We  will  not  consider  m  or  Yq  as  they  occur  only  in  the  numerator  of  Equation  (2.28) 
and  obviously  maximize  power  by  taking  on  the  largest  value  possible: 


P  C(g)3^3 

mio2  (l-g)2  +  (2Cf)2 


(4.1) 


The  second  model  is  an  adaptation  of  Equation  (3.39).  It  consists  of  six  variables.  We 
will  not  consider  m,  y,  or  H  as  they  occur  only  in  the  numerator  of  Equation  (3.39)  and 
obviously  maximize  power  by  taking  on  the  largest  value  possible: 

P  _  _ A2RCObaSe6{  1  +Cq2 CObase2R2)2 _ 

( myH )  2[((—  (0/,rtAY.-  r  I  \-Cq  CObase2R2')-\  CoA2(Obase2R2')2  \-  (2wb(_;  ^  I  ^-Cq  Q)base~R~)-\-A2G)baseR')2 

(4.2) 

Equations  (4.1)  and  (4.2)  lie  in  the  core  of  our  analysis. 


35 


4.3  Variable  Ranges 

In  order  to  optimize  maximum  power  generated  we  must  constrain  our  variables.  The  goal 
of  the  optimization  is  to  find  the  optimum  combination  of  variables  in  order  to  maximize 
power  and  lead  to  the  development  of  more  efficient  piezoelectric  generators.  Keeping  in 
mind  the  practical  application  of  our  optimization,  we  select  the  ranges  of  each  variable  in 
such  a  way  that  engineers  will  be  able  to  select  materials  and  conditions  within  the  limits  of 
each  variable  to  create  a  physical  model  based  on  our  results.  Figures  4.1  and  4.2  display 
the  ranges  for  variables  in  Equations  (4.1)  and  (4.2),  respectively. 

Table  4.1:  Variable  ranges  for  the  three-variable  power  generation  model  of  Equation  (4.1). 


Variable 

Variable  Symbol 

Range 

Frequency  of  Excitation 

ft) 

[120^“1,360^-1] 

Natural  Frequency 

(On 

[120^"1,360^-1] 

Modal  Damping  Ratio 

c 

[0.005,0.02] 

In  Roundy  and  Wright’s  experiment  [15],  they  use  the  values  ft)  =  (On  =  240 ns~l  and  £  = 
0.015.  We  choose  the  range  of  ft) ,  ft),,  and  £  by  varying  Roundy  and  Wright’s  experimental 
values. 


Table  4.2:  Variable  ranges  for  the  six-variable  power  generation  model  of  Equation  (4.2). 


Variable 

Variable  Symbol 

Range 

Electromechanical  Coupling  Coefficient 

A 

[0.01,0.99] 

Load  Resistance 

R 

[0.02H,  100f2] 

Base  Motion  Driving  Frequency 

(&base 

[1207T^“1,3607r^-1] 

Net  Clamped  Cap.  of  Piezoelectric  Material 

Co 

[1  x  10”8F,  1  x  10”6F] 

Modal  Short  Circuit  Net  Frequency 

Wsc 

[12071V”1, 36071V”1] 

Modal  Damping  Ratio 

c 

[0.005,0.02] 

We  choose  the  range  of  CObase  and  £  by  varying  the  experimental  values  found  in  Roundy 
and  Wright’s  paper  [15].  Likewise,  we  choose  the  range  of  Co  by  varying  the  value  cited 
in  Fleming  and  Behren’s  paper  [24].  We  choose  to  use  the  same  range  for  wsc  as  the 
other  frequency,  co^ase ■  For  the  coefficient,  A,  we  use  [0.01,0.99].  Finally,  the  range  of 
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R  varies  from  the  resistance  of  a  copper  wire  to  about  half  the  resistance  of  an  operating 
incandescent  light  bulb. 


4.4  mk  Factorial  Sampling 

One  method  to  find  an  optimal  solution  for  power  generation  using  our  models  is  to  evaluate 
the  combination  of  all  independent  variables,  which  we  will  interchangeably  refer  to  as 
factors  in  this  chapter,  at  different  levels.  Levels  are  the  variety  of  values  a  factor  may 
assume  in  a  simulation.  This  methodology  can  be  described  as  mk  factorial  design,  where 
m  is  the  number  of  levels  per  factor  and  k  is  the  number  of  factors.  In  an  nr  factorial 
design  the  number  of  design  points  is  equal  to  mk .  Design  points  are  the  combination  of 
factor  levels  to  be  evaluated  in  the  simulation.  The  larger  the  value  of  m  for  an  nr  factorial 
design,  the  better  its  space-filling  properties  [25].  However,  an  increase  of  levels  may 
drastically  increase  the  number  of  design  points  required  to  evaluate  in  a  simulation. 


4.5  Nearly  Orthogonal  and  Balanced  Latin  Hypercube 
Sampling 

An  alternate  and  more  efficient  sampling  methodology  is  to  use  space-filling  designs  to 
thoughtfully  choose  the  specific  levels  of  each  factor  to  evaluate.  Latin  hypercube  designs 
provide  a  flexible  way  of  constructing  efficient  designs  for  qualitative  factors.  They  have 
some  of  the  space-filling  properties  of  mk  factorial  design,  but  require  orders  of  magnitude 
less  sampling  [25]. 

K.Q.  Ye  [1998]  describes  Latin  hypercubes  as  follows: 

An  n  x  d  Latin  hypercube  consists  of  d  permutations  of  Sn  =  {1,2 
A  Latin  hypercube  design  takes  row  vectors  as  the  experimental  points  in 
a  J-dimensional  space.  One-dimensional  projections  of  a  Latin  hypercube 
design  are  evenly  spaced  and  have  the  lowest  possible  discrepancy.  [26] 

The  significant  benefits  of  the  Latin  hypercube  sampling  used  in  this  paper  is  the  samplings’ 
orthogonality  and  space-filling  property. 
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4.6  Sampling  Metric:  Orthogonality 

The  correlation  between  two  vectors  u  —  [u\,u.2,  ...un\  and  v  =  [vi,V2,...v„]  is  defined  the 
value  of  the  quotient 

L(vi-v)(uj-u)  (43) 

a/L(m/-m)2L(v;--v)2’ 

where  u  —  ^  and  v  —  If  u  and  v  are  uncorrelated  they  are  said  to  be  orthogonal 
and  their  correlation  value  is  zero.  Alternatively,  two  completely  correlated  vectors  have 
a  correlation  value  of  one.  Orthogonal  Latin  hypercubes  create  column  vectors  with  zero 
correlation  [26]. 

Orthogonality  in  variable  sampling  is  particularly  valuable  in  regression  analysis  and  sim¬ 
ulations  where  one  desires  independent  variables  in  a  regression  model  to  be  orthogonal  to 
each  other  in  order  to  minimize  variable  correlation  [27],  [28].  Additionally,  the  orthog¬ 
onal  Latin  hypercube  design  ensures  independence  of  estimates  of  linear  effects  of  each 
variable  [26]. 

No  matter  how  many  levels  or  factors,  nr  factorial  samplings  result  in  completely  orthog¬ 
onal  design  points.  In  other  words,  the  design  columns  in  factorial  design  are  uncorrelated. 

We  use  a  select  set  of  columns  from  a  nearly  orthogonal  and  balanced  (NOB)  design, 
developed  by  H.  Vieira  Jr.  [29],  in  order  to  reduce  the  number  of  design  points  necessary  to 
efficiently  find  the  optimal  solution  to  our  power  generation  problem.  After  we  input  our 
model  variables,  with  associated  variable  ranges,  Vieira’s  NOB  design  spreadsheet  outputs 
512  nearly  orthogonal  and  balanced  design  points  [30].  According  to  Vieira, 

Nearly  orthogonal  means  that  the  maximum  absolute  pairwise  correlation 
between  and  two  design  columns  is  minimal.  Nearly  balanced  means  that 
for  any  single  factor  column,  the  number  of  occurrences  of  each  distinct 
factor  level  is  nearly  equal. 

Vieira  and  Sanchez  call  a  design  nearly  orthogonal  if  the  absolute  pairwise  correlation 
between  any  two  factors  is  less  than  0.05.  Likewise,  they  call  a  design  nearly  balanced  if 
the  number  of  objects  in  each  level  of  each  factor  differs  from  the  ideal  by  no  more  than 
20%  [25], 
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The  balance  is  a  concern  when  discrete-valued  factors  with  limited  numbers  of  levels  are 
present  because  rounding  the  levels  appropriate  for  continuous-valued  factors  may  induce 
correlations  that  are  larger  than  desired.  Because  the  designs  we  use  in  this  paper  involve 
only  continuous-valued  factors,  hence  the  subset  of  columns  is  a  Latin  hypercube,  we  will 
refer  to  these  as  nearly  orthogonal  and  balanced  Latin  hypercube  (NOBLH)  designs. 

Below  are  the  correlation  between  factors  for  the  nearly  orthogonal  and  balanced  Latin 
hypercube  sampling,  consisting  of  the  5 12  design  points  generated  by  Vieira’s  NOB  Design 
Spreadsheet.  Refer  to  Appendix,  Section  A.l  and  Section  A. 2  to  see  the  MATLAB  code 
we  used  to  generate  orthogonality  values  cited  in  Tables  4.3  and  4.4. 


Table  4.3:  Orthogonality  of  three-variable  NOB  Latin  hypercube  design. 


ft) 

C 

ft) 

-0.0011 

-0.005 

-0.0011 

-0.0038 

c 

-0.005 

-0.0038 

Table  4.4:  Orthogonality  of  six-variable  NOB  Latin  hypercube  design. 


A 

R 

®base 

Co 

Wsc 

c 

A 

-0.00106 

-0.0036 

-0.00293 

0.001528 

0.001367 

R 

-0.00106 

-0.0038 

-0.00267 

-0.00195 

-0.00169 

®base 

-0.00359 

-0.00378 

-0.0017 

0.002429 

0.001984 

Co 

-0.00293 

-0.00267 

-0.0017 

-0.00253 

-0.00179 

Wsc 

0.001528 

-0.00195 

0.00243 

-0.00253 

0.006878 

c 

0.001367 

-0.00169 

0.00198 

-0.00179 

0.006878 

Recall  that  perfect  orthogonality  is  defined  as  having  a  correlation  value  of  zero.  Notice  the 
minimum  absolute  pairwise  correlation  between  each  design  column.  Using  our  variable 
inputs,  Vieira’s  NOBLH  design  produces  close  to  orthogonal  design  columns  for  both  our 
three-  and  six- variable  models. 
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4.7  Sampling  Metric:  Space-Filling  Property 


It  is  desirable  to  select  sample  design  points  with  good  space-filling  properties,  in  other 
words,  design  points  distributed  throughout  the  entire  experimental  region.  For  this  paper 
we  rate  the  space-filling  of  the  samplings  using  the  Euclidean  Maximin  (Mm)  distance. 

For  a  given  design,  define  a  distance  list  d  =  (d\,d.2,  ...,drn(n_i)i/2),  where 
the  elements  of  d  are  the  Euclidean  inter-site  distance  of  the  n  design  points, 
ordered  from  smallest  to  largest.  The  Euclidean  Mm  distance  is  defined  as 
di,  where  a  larger  value  is  better.  A  large  value  of  d\  means  that  no  two 
points  are  close  to  (within  d\  of)  each  other.  [27] 

We  use  the  Equation  (4.4)  to  compute  the  Euclidean  Mm  distance: 

k 

d(s,  t)=  Y,  sj~TJ 2  (4-4) 

V-1 

where  the  design  space  is  an  n  x  k  matrix  and  s  and  t  are  any  two  design  points  [31]. 

When  comparing  the  Euclidean  Mm  distance  of  two  samples  with  equal  number  of  design 
points,  it  is  possible  to  depend  solely  on  Equation  (4.4).  However,  nr  factorial  and  NOBLH 
designs  result  in  differing  numbers  of  design  points.  NOBLH  designs  result  in  512  design 
points,  while  the  number  of  mk  factorial  design  points  is  a  function  of  levels  (m)  and  factors 
(k).  For  this  Euclidean  Mm  distance  calculation  we  use  a  three -variable  factorial  design 
with  10  levels  and  a  six-variable  factorial  design  with  five  levels,  resulting  in  1000  and 
15,625  design  points,  respectively. 

In  order  to  compare  Euclidean  Mm  distance  of  NOBLH  and  mk  factorial  designs,  we  first 
normalize  the  design  columns  of  each  design  with  values  from  the  mk  factorial  design. 
After  column  normalization,  it  is  possible  to  use  Equation  (4.4)  to  compute  the  Euclidean 
Mm  distance  and  compare  the  results.  Refer  to  Appendix,  Section  A. 3  for  our  Euclidean 
Maximin  Distance  Code. 

The  three-variable  and  six-variable  NOBLH  designs  result  in  Euclidean  Mm  distances  of 
0.018  and  0.1947,  respectively.  The  three-  and  six- variable  mk  factorial  designs  result 
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in  Euclidean  Mm  distances  of  0.1111  and  0.25,  respectively.  In  both  the  three-  and  six- 
variable  models,  the  mk  factorial  designs  provide  space-filling  properties  optimal  to  the 
NOBLH  designs. 

Of  course,  another  way  of  judging  the  space-filling  property  of  a  design  is  by  observing 
the  location  of  the  design  points  in  reference  to  the  factor  ranges.  Unfortunately,  this  vi¬ 
sual  technique  is  only  possible  in  three  dimensions  or  below.  See  Figure  4. 1  for  a  visual 
depiction  of  the  1000  design  points  evaluated  using  10  levels  (in)  for  the  three-factor  (k) 
model. 


1100 


Figure  4.1:  mk  factorial  sampling. 
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Figure  4.2  visually  depicts  the  512  nearly  orthogonal  and  balanced  Latin  hypercube  design 
points  evaluated  for  the  three- variable  model. 


Figure  4.2:  Nearly  orthogonal  and  balanced  Latin  hypercube  sampling. 


As  the  reader  will  notice,  the  NOBLH  sampling  method  samples  extremely  well  from  the 
interior  of  the  design  space  and  does  a  poor  job  sampling  from  the  exterior  of  the  design 
space.  The  mk  factorial  sampling  method,  on  the  other  hand,  samples  uniformly  across  the 
design  space. 
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4.8  Sampling  Metric:  Efficiency 

After  evaluating  all  design  points  in  our  model,  we  determine  the  maximum  power  gener¬ 
ated  and  determine  an  optimal  combination  of  factor  levels.  This  nr  factorial  method  to 
evaluate  the  six- variable  model  is  easily  accomplished  using  the  MATLAB  code  referenced 
in  Appendix,  Section  A. 6. 

The  six-variable  power  generation  model  requires  49  floating  point  operations  (flops)  in 
order  to  compute  one  design  point.  The  MATLAB  code  in  Appendix,  Section  A. 6  evaluates 
the  power  generated  for  all  combinations  of  the  six  factors  at  a  varying  number  of  levels. 
This  equates  to  m6  design  points.  Therefore  the  efficiency  of  the  algorithm  with  the  given 
number  of  factors  and  levels  is  49  x  m6  flops.  This  model  only  consists  of  six  factors  and 
therefore  the  number  of  design  points  required  does  not  change  as  dramatically  with  an 
increase  of  levels  as  a  model  with  more  factors.  For  this  paper,  we  use  as  many  as  15  levels 
per  factor.  When  evaluating  the  model  with  15  levels  per  factor,  156  =  11, 390, 625  design 
points,  and  49  x  1 1, 390, 625  ~  558  million  flops,  our  Core  i7  equipped  computer  can  carry 
out  the  algorithm  in  just  over  two  minutes. 

If  our  model  consisted  of  ten  factors,  an  nr  factorial  sampling  of  15  levels  would  consist 
of  1510^  5.77  x  1011  design  points  and  too  many  flops  for  a  standard  desk  top  computer 
to  execute  in  a  week.  This  algorithm  and  mk  factorial  design  are  not  a  wise  use  of  time  for 
simulations  with  a  large  number  of  factors,  especially  if  we  desire  to  change  the  any  aspect 
of  the  factors  or  levels  and  reproduce  the  simulation. 

The  real  benefit  of  nearly  orthogonal  and  balanced  Latin  hypercube  design  is  its  efficiency. 
Using  Vieira’s  NOB  design  it  is  possible  to  evaluate  a  model  of  up  to  20  k-level  discrete 
factors  (where  k  =  2, 3...  1 1  levels)  and  100  continuous  factors  using  only  512  design  points 
[30].  Many  models  or  simulations  contain  many  more  factors  and  computations.  NOB 
designs  choose  very  good  samplings  to  test  given  the  limited  number  of  simulations  a 
system  can  compute  in  a  given  time. 

4.9  NOBLH  Iterative  Method 

An  additional  optimization  method  involves  multiple  iterations  of  sampling  the  design 
space  and  using  results  to  restrict  the  design  space  to  choose  design  points.  In  this  pa- 
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per  we  refer  to  this  method  as  the  NOBLH  Iterative  method. 

First,  we  use  NOBLH  design  to  choose  design  points  from  the  original  design  space  set  by 
our  variable  ranges.  Next,  we  analyze  the  design  point  that  generate  the  maximum  power 
output.  By  studying  the  values  of  each  variable  in  the  optimal  design  point,  we  are  able 
to  make  educated  guesses  about  how  to  restrict  the  variable  ranges  in  order  to  focus  our 
design  space. 

With  each  newly  discovered  optimal  design  point  we  are  able  to  restrict  the  variable  ranges 
and  design  space  until  we  reach  a  maximum  power  or  identify  an  issue  that  causes  a  cease 
to  iterations.  The  aspiration  of  this  method  is  to  arrive  at  similar  results  to  that  of  the  more 
costly  mk  factorial  sampling  method  while  evaluating  fewer  design  points. 

It  should  be  noted,  the  iterative  method  we  use  is  not  as  systematic  as  other  iterative  meth¬ 
ods.  Through  the  use  of  our  NOBLH  iterative  method  we  wish  to  gain  insights  into  our 
optimization  models  and  better  understand  alternative  methods  to  optimization. 
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CHAPTER  5: 

Optimization  Results  and  Analysis 


5.1  Three-Variable  Model  NOBLH  and  mk  Factorial 
Results 

After  sampling  our  design  space  using  both  NOBLH  and  nr  factorial  sampling  methods  we 
use  the  resulting  design  points  to  find  the  optimal  values  of  the  three  variables  (ft),  ft),, ,  and 
£)  in  order  to  maximize  power  generation.  For  the  nr  factorial  samplings  we  use  levels 
of  10  and  100.  See  the  Appendix  for  the  Three- Variable  Model  Optimization  Codes  for 
factorial  sampling  (Section  A.4)  and  NOBLH  sampling  (Section  A. 5).  The  results  of  each 
sampling  method  are  summarized  in  Table  5.1. 


Table  5.1:  Optimization  of  power  in  the  three-variable  model. 


Sampling  Method 

Design  Points 

Optimal  Power 

ft) 

ft),, 

c 

NOBLH 

512 

2.49  x  1010 

339.87T 

339.87T 

0.0122 

mk  (m  =10) 

1,000 

7.23  x  1010 

3607T 

360;r 

0.005 

mk  (, m  -  100) 

1,000,000 

7.23  x  1010 

360; r 

360;r 

0.005 

The  reader  will  notice  both  nf  factorial  sampling  methods  are  able  to  find  an  optimal 
power  close  to  three  times  larger  than  the  NOBLH  sampling  method.  The  reader  will  also 
notice  both  nr  factorial  samplings  determine  the  optimal  values  of  each  factor  to  be  at  the 
extreme  of  the  factors’  ranges.  The  optimal  values  of  ft)  and  ft),,  are  the  maximum  values 
allowed  by  the  range,  3607T.  and  the  optimal  value  of  £  is  the  minimum  value  allowed  by 
the  range,  0.005.  As  noted  in  Chapter  4,  the  NOBLH  sampling  fails  to  adequately  sample 
the  extremes  of  the  design  space,  while  the  mk  factorial  method  samples  the  entire  design 
space  uniformly. 

This  failure  of  sampling  design  points  in  the  extremities  of  the  design  space  accounts  for 
only  part  of  the  discrepancy  of  optimal  solutions  discovered  by  the  sampling  techniques. 
Notice  ft)  =  ft),,  for  both  optimal  values  for  the  nr  factorial  samplings  and  the  NOBLH 
sampling.  After  further  investigation  into  the  model,  we  find  the  optimal  values  for  the  ft) 
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and  con  are  always  equal. 


We  can  prove  this  mathematically  by  analyzing  Equation  (4.1).  In  order  to  maximize  the 
model  we  obviously  want  to  make  the  numerator  as  large  as  possible  while  making  the 
denominator  as  small  as  possible.  In  order  to  eliminate  the  first  part  of  the  denominator, 
(1  —  ^t)2,  we  set  co  =  con.  This  reduces  the  model  to  the  following  equation: 


p  c  ®3  «3 

mY^  ~  4^"  4£' 


(5.1) 


The  three-variable  model  reduces  to  a  two-variable  model  in  which  co  —  con.  In  order  to 
maximize  power,  one  simply  maximizes  co  and  minimizes  £. 


After  analyzing  the  design  points,  we  discover  that,  out  of  the  512  design  points  in  the 
NOBLH  sampling,  only  four  design  points  set  co  =  C0n,  compared  to  100  and  1000  points 
for  the  m  =  10  and  m  =  100  nr  factorial  samplings,  respectively. 
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Setting  co  —  (On  plays  a  crucial  role  in  optimizing  this  model.  The  m  =10  mk  factorial 
design  selects  100  design  points  with  co  —  (On.  Out  of  the  said  100  design  points,  16  de¬ 
sign  points  produce  power  greater  than  the  maximum  power  generated  by  any  of  the  512 
NOBLH  design  points.  See  Figure  5.1  for  the  sixteen  mr  factorial  design  points  that  out¬ 
perform  the  NOBLH  sampling. 


Figure  5.1:  mk  factorial  design  points  achieving  higher  power  generation  than  all  NOBLH  design 
points. 


While  optimizing  this  particular  model,  mr  factorial  sampling  proves  the  best  sampling 
method.  After  analysis,  we  determine  that  to  maximize  power  one  can  set  co  =  (On  and 
minimize 
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5.2  Three- Variable  Model  NOBLH  Iterative  Results 

After  determining  the  optimal  power  generation  and  associated  variable  values  for  the 
three-variable  model  through  the  mr  factorial  results,  we  attempt  to  replicate  the  results 
using  the  NOBLH  iterative  method.  Our  initial  hopes  are  the  NOBLH  iterative  method 
would  be  able  to  reach  similar  results  while  evaluating  less  design  points  than  the  more 
costly  nr  factorial  sampling  method.  In  this  particular  case,  the  three  level  nr  factorial 
sampling  finds  the  optimal  power  using  only  1000  design  points.  Even  if  the  NOBLH  iter¬ 
ative  method  finds  the  optimal  results  in  two  iterations,  the  method  would  use 
512  x  2  =  1, 024  design  points  to  do  so. 

The  first  iteration  of  the  NOBLH  iterative  method  is  identical  to  the  NOBLH  sampling 
method.  Table  5.2  uses  the  same  variable  ranges  as  in  Table  4.1  and  Table  5.3  arrives  at  the 
same  results  as  in  Table  5.1. 

Table  5.2:  Iteration  one  NOBLH  iterative  variable  ranges. 


ft) 

ft)/1 

c 

Lower  Limit 

1207T 

1207T 

0.005 

Upper  Limit 

3607T 

3607T 

0.02 

Table  5.3:  Iteration  one  NO 

BLH  iterative  results. 

Optimal  Power 

ft) 

(On 

C 

2.49  x  1010 

339.Sk 

339.8  K 

0.01222 

After  analyzing  the  results  listed  in  Table  5.3,  we  are  able  to  restrain  the  ranges  of  ft),  (On 
and  £,  listed  in  Table  5.4.  Using  the  new  variable  ranges  we  produce  512  new  design 
points  using  NOBLH  sampling.  After  evaluating  the  new  design  points,  we  determine  the 
resulting  optimal  power  generation  and  associated  variable  values,  as  shown  in  Table  5.5. 


Table  5.4:  Iteration  two  NOBLH  iterative  variable  ranges. 


ft) 

(On 

c 

Lower  Limit 

3 18.3zr 

318.3  K 

0.005 

Upper  Limit 

360 n 

360 K 

0.015 
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Table  5.5:  Iteration  two  NOBLH  iterative  results. 


Optimal  Power 

ft) 

(On 

C 

5.38  x  1010 

349.27T 

349.97T 

0.0055 

Analyzing  the  results  in  Table  5.5,  we  see  optimal  power  increase  again  along  with  both 
values  of  ft)  and  ft),,.  The  optimal  value  of  £  continues  to  decrease.  We  use  this  data  to 
further  tighten  all  variable  ranges  for  the  third  iteration,  listed  in  Table  5.6. 

Table  5.6:  Iteration  three  NOBLH  iterative  variable  ranges. 


ft) 

ft),, 

c 

Lower  Limit 

340.67 t 

340.67T 

0.005 

Upper  Limit 

3607 t 

3607 t 

0.006 

Table  5.7:  Iteration  three  NOBLH  iterative  results. 


Optimal  Power 

ft) 

(On 

C 

6.73  x  1010 

357.  Itt 

357  An 

0.005135 

Table  5.7  displays  an  additional  increase  in  optimal  power,  ft),  and  ft),,.  Additionally,  the 
optimal  value  of  £  continues  to  decrease.  We  are  able  to  restrain  the  ranges  for  all  variables 
for  the  fourth  iteration,  as  listed  in  Table  5.8. 

Table  5.8:  Iteration  four  NOBLH  iterative  variable  ranges. 


ft) 

(On 

c 

Lower  Limit 

353.37T 

353.37T 

0.005 

Upper  Limit 

3607T 

3607 t 

0.0052 

Table  5.9:  Iteration  four  NOBLH  iterative  results. 


Optimal  Power 

ft) 

(On 

C 

7.12  x  1010 

359.47T 

359.57T 

0.005029 

Table  5.9  displays  optimal  power  close  to  the  results  found  using  mk  factorial  sampling. 
After  analysis,  we  further  restrain  variable  ranges  for  the  fifth  iteration,  listed  in  Table 
5.10. 
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Table  5.10:  Iteration  five  NOBLH  iterative  variable  ranges. 


ft) 

(On 

c 

Lower  Limit 

359; t 

359k 

0.005 

Upper  Limit 

360;r 

360 n 

0.0052 

able  5.11:  Iteration  five  NOBLH  iterative  results. 


Optimal  Power 

ft) 

(On 

c 

7.23  x  1010 

369k 

360 n 

0.005002 

After  five  iterations  we  find  optimal  power  and  associated  variable  values  (listed  in  Table 
5.11)  equivalent  to  mk  factorial  results.  Our  conservative  methodology  to  iterative  range 
restriction  forces  us  to  reach  our  results  after  analyzing  5x512  =  2, 560  design  points.  The 
three  level  mk  factorial  sampling  optimization  method  determines  the  same  results  after 
computing  only  1000  design  points.  However,  the  NOBLH  iterative  method  shows  great 
potential  when  dealing  with  a  large  number  of  factors. 


5.3  Six- Variable  Model  NOBLH  and  mk  Factorial 

Results 

After  sampling  our  design  space  using  both  NOBLH  and  nr  factorial  sampling  methods, 
we  use  the  resulting  design  points  to  find  the  optimal  values  of  the  six  variables  (A,  R, 
( Obase >  Cq,  wsc  and  £)  in  order  to  maximize  power  generation.  For  mk  factorial  samplings 
we  use  levels  of  3,  5,  7,  10,  11,  and  12.  See  the  Appendix  for  the  Six- Variable  Model 
Optimization  Codes  for  factorial  sampling  (Section  A. 6)  and  NOBLH  sampling  (Section 
A.7).  The  results  of  each  sampling  method  are  given  in  Table  5.12. 
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Table  5.12:  Optimization  of  power  in  the  six-variable  model. 


Method 

Design  Points 

Optimal  Power 

A 

R 

®base 

Co 

Wsc 

£ 

NOBLH 

512 

7.8509  x  109 

0.754 

13.129 

342.61tt 

4.1  x  10~8 

344.96 n 

0.006 

mk  (m  =  3) 

729 

1.8082  x  1010 

0.5 

50.01 

360;r 

1  x  10-6 

360 k 

0.005 

mk  (m  =  5) 

15,625 

1.8082  x  1010 

0.5 

50.01 

360tz; 

1  x  10~6 

360 n 

0.005 

mk  (m  =  7) 

117,649 

1.8256  x  1010 

0.337 

100 

360tt 

1  x  10-6 

360 n 

0.005 

mk  ( m  =10) 

1,000,000 

1.8256  x  1010 

0.337 

100 

360 k 

1  x  10~6 

360 n 

0.005 

mk  (m  =  11) 

1,771,561 

1.8168  x  1010 

0.402 

70.01 

360 k 

1  x  10~6 

360 k 

0.005 

mk  (m  =  12) 

2,985,984 

1.8203  x  1010 

0.356 

90.91 

360 n 

1  x  10~6 

360 k 

0.005 

Like  the  three-variable  model,  all  nr  factorial  samplings  discover  combinations  of  factor 
values  capable  of  generating  power  more  than  twice  the  value  of  the  NOBLH  sampling’s 
optimal  solution.  Also  like  the  three-variable  model,  the  optimal  solution  in  the  mk  factorial 
samplings  includes  four  variables  at  the  extremes  of  their  ranges.  The  optimal  values  of 
Cdbase,  Qb  and  wsc  are  the  maximums  of  the  their  ranges,  while  the  optimal  value  for  £  is 
the  minimum  of  its  range.  As  we  noted  in  the  previous  section,  NOBLH  sampling  fails  to 
sample  the  extremes  of  the  design  space.  The  simulations  using  the  NOBLH  design  points 
does  not  have  the  opportunity  to  explore  power  generation  at  the  extremes  of  the  design 
space.  With  only  217  more  design  points  than  the  NOBLH  sampling,  the  three  level  nr 
factorial  sampling  creates  design  points  in  the  extremes  of  the  variable  ranges  and  finds  an 
optimal  solution  more  than  twice  the  value  of  the  NOBLH  design  solution. 

An  interesting  result  of  our  simulation  is  that  creating  a  finer  grid  (increasing  the  number 
of  levels  and  resulting  design  points)  does  not  always  result  in  a  more  optimal  solution. 
For  example,  both  the  m  =  7  and  m  =10  mk  factorial  samplings  find  an  optimal  solution 
of  1.8256  x  1010.  However,  when  the  number  of  levels  increases  from  ten  to  eleven  the 
optimal  solution  decreases  to  1.82168  x  1010.  Additionally,  the  reader  will  notice  the  max¬ 
imum  power  found  with  m  =12  sampling  is  also  less  than  that  discovered  with  m  =10 
sampling.  This  can  be  explained  by  the  fact  both  the  m  =  7  and  m  =  10  samplings  include 
A  =  0.337,  while  the  m  =11  and  m  =  12  samplings  do  not. 

The  true  optimal  value  depends  on  the  precise  sampling  of  A  and  R.  Using  the  assumption 
the  optimal  power  generation  value  occurs  at  the  maximum  allowable  values  of  (Obase ,  Qb 
and  wsc  and  minimum  allowable  value  of  £,  we  are  able  to  further  refine  the  true  optimal 
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values  of  A  and  R.  We  accomplish  this  refinement  by  fixing  (Obase,  Cq,  wsc,  and  £  at  their 
optimal  value  while  sampling  1000  levels  each  of  A  and  R.  This  results  in  refining  the 
optimal  values  of  A  as  0.3376  and  R  as  100.  With  this  more  precise  value  for  A  we  arrive  at 
a  new  optimal  power  generation  only  6.73  x  104  greater  than  our  previous  optimal  solution. 

5.4  Six- Variable  Model  NOBLH  Iterative  Results 

After  determining  the  optimal  power  generation  and  associated  variable  values  for  the  six- 
variable  model  through  the  mr  factorial  results,  we  attempt  to  replicate  the  results  using  the 
NOBLH  iterative  method.  Our  aspiration  is  the  NOBLH  iterative  method  will  be  able  to 
reach  similar  results  while  evaluating  less  design  points  than  the  more  costly  mk  factorial 
sampling  method.  From  Table  5.12,  we  see  using  seven  levels  the  nr  factorial  sampling 
method  results  with  an  optimal  power  of  1.8256  x  1010  after  evaluating  117,649  design 
points.  We  deem  the  NOBLH  iterative  method  a  success  if  we  reach  a  optimal  power 
close  to  the  seven  level  mr  factorial  sampling  method  results  with  far  fewer  design  points 
evaluated. 

The  first  iteration  of  the  NOBLH  iterative  method  is  identical  to  the  NOBLH  sampling 
method.  Table  5.13  uses  the  same  variable  ranges  as  in  Table  4.2,  and  Table  5.14  arrives  at 
the  same  results  as  in  Table  5.12. 


Table  5.13:  Iteration  one  NOBLH  iterative  variable  ranges. 


A 

R 

®base 

Co 

Wsc 

c 

Lower  Limit 

0.01 

0.02 

120  n 

1  x  10~8 

120; r 

0.005 

Upper  Limit 

0.99 

100 

360;r 

1  x  10”6 

360;r 

0.02 

Table  5.14 

:  Iteration 

one  NOBI 

LH  iterative  results. 

Optimal  Power 

A 

R 

®base 

Co 

wsc 

C 

7.85  x  109 

0.754 

13.129 

342. 6;r 

4.1  x  10~8 

345k 

0.005998 

After  analyzing  the  results  listed  in  Table  5.14,  we  are  able  to  assume  optimal  power  in¬ 
creases  when  A,  R,  (Obase  and  w'sc  take  on  values  closer  to  their  upper  most  range  value. 
We  also  assume  £  and  Co  taking  values  closer  to  their  minimum  allowable  values  results  in 
an  increase  in  power.  From  these  assumptions  we  tighten  the  ranges  of  all  three  variables, 
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listed  in  Table  5.15.  Using  the  new  variable  ranges  we  produce  512  new  design  points  using 
NOBLH  design.  After  evaluating  the  new  design  points  we  determine  the  resulting  optimal 
power  generation  and  associated  variable  values,  given  in  Table  5.16. 


Table  5.15:  Iteration  two  NOBLH  iterative  variable  ranges. 


A 

R 

®base 

Co 

Wsc 

c 

Lower  Limit 

0.25 

1 

340k 

1  x  10~8 

340k 

0.005 

Upper  Limit 

0.99 

100 

360;r 

1  x  10~7 

360 n 

0.01 

Table  5.16:  Iteration  two  NOBLH  iterative  results. 


Optimal  Power 

A 

R 

t &base 

Co 

WSC 

C 

1.64  x  1010 

0.812 

13.98 

358.57T 

1.28  x  10-8 

358.7tt 

0.00533 

After  the  second  iteration  we  notice  an  increase  in  optimal  power,  (0base,  and  wsc-  The 
optimal  A  and  R  increase  slightly.  The  optimal  £  and  Co  reduce  slightly.  Using  these 
findings  we  further  restrain  variable  ranges,  listed  in  Table  5.17.  Our  results  are  shown  in 
Table  5.18. 


Table  5.17:  1 

teration  three  NOBL 

H  iterative  variable  ranges. 

A 

R 

t fibase 

Co 

wsc: 

c 

Lower  Limit 

0.5 

5 

3587T 

1  x  10~8 

3587T 

0.005 

Upper  Limit 

0.99 

50 

360; r 

5  x  10~8 

360k 

0.006 

Ta 

ble  5.18: 

Iteration  three  NOBLH  iterative 

results. 

Optimal  Power 

A 

R 

®base 

Co 

wsc 

c 

1.778  x  1010 

0.626 

28.24 

358. 4;r 

4.1  x  10~8 

35%.2k 

0.005006 

From  the  third  iteration  we  notice  an  increased  optimal  power  along  with  a  slight  decrease 
in  the  optimal  values  of  (Obase,  wsc,  and  £.  The  optimal  values  of  A,  R ,  and  Co  increase 
slightly.  We  continue  restricting  variable  ranges,  listed  in  Table  5.19. 
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Table  5.19: 

Iteration  four  NOB 

LH  iterative 

variable 

ranges. 

A 

R 

t &base 

Co 

wsc 

c 

Lower  Limit 

0.5 

20 

3587T 

2  x  10'8 

358tt 

0.005 

Upper  Limit 

0.8 

30 

359; t 

5  x  10”8 

359k 

0.00501 

Tab 

le  5.20:  1 

Iteration  four  NOBI 

LH  iterative  results. 

Optimal  Power 

A 

R 

(fibase 

Co 

Wsc 

c 

1.792  x  1010 

0.728 

21.31 

358. 9;r 

2.9  x  10~8 

358. 9;r 

0.005 

From  the  results  of  the  fourth  iteration,  Table  5.20,  we  see  another  increase  in  optimal 
power.  However,  knowing  the  optimal  values  discovered  from  nr  factorial  sampling,  we 
have  made  incorrect  assumptions  restricting  the  ranges  of  A,  R,  (Obase,  Co,  and  wsc-  If 
we  further  restrict  the  variable  ranges  in  subsequent  iterations  we  would  fail  to  reach  an 
optimal  power  similar  to  that  found  through  nr  factorial  sampling. 

An  alternative  way  to  look  at  the  NOBLH  iterative  results  is  purely  in  terms  of  the  number 
of  design  points  evaluated.  After  four  iterations,  we  only  evaluate  4x512  =  2048  design 
points.  The  optimal  power  resulting  from  the  NOBLH  iterative  method  (1.792  x  1010) 
only  represents  a  1.84%  loss  in  power  from  the  seven  level  nr  factorial  sampling  method 
(1.8256  x  1010).  Less  than  a  two  percent  loss  in  optimal  power  seems  minor  when  you 
consider  the  seven  level  mk  factorial  sampling  method  requires  the  evaluation  of  117,649 
design  points.  Again,  the  NOBLH  iterative  method  shows  potential  when  optimizing  mod¬ 
els  with  a  large  number  of  variables. 

5.5  Robustness  Analysis 

The  sampling  methods  provides  the  means  to  find  the  optimal  combination  of  variable 
values  to  maximize  power  generation  in  our  models.  After  running  simulations  with  the 
sampled  design  points  we  find  the  optimal  solution  for  each  model.  On  paper,  we  then 
possess  the  exact  metrics  to  build  a  piezoelectric  generator  prototype  in  order  to  provide 
optimal  power  output.  However,  the  precision  of  the  variable  values  we  insist  on  paper 
are  not  always  possible  to  replicate  in  experiments  or  practical  application,  even  under 
perfect  conditions.  All  variables  depending  on  materials  or  conditions  have  associated 
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variability.  For  example,  a  resistor  may  be  advertised  as  providing  a  certain  number  of 
ohms  of  resistance,  but  in  reality,  provides  the  advertised  ohms  of  resistance  plus  or  minus 
some  level  of  discrepancy  allowed  by  the  manufacturer. 

The  purpose  of  robust  design,  as  a  process  of  simulation  optimization,  is  to  identify  the 
“best”  solution  as  one  not  overly  sensitive  to  small  changes  in  system  inputs.  In  this  section 
we  will  attempt  to  identify  which  factors  in  the  design  of  a  piezoelectric  generator  need 
more  precision  and  which  factors  will  allow  more  variability  while  producing  as  much 
power  as  possible  [32]. 

5.5.1  Three-Variable  Model 

Using  mathematical  models  to  simulate  a  real  world  problem,  it  is  easy  to  see  the  variables 
only  as  their  numerical  values.  In  reality  the  variables  respond  to  specifications  of  the 
piezoelectric  generator  and  its  environmental  setting. 

In  this  simple  three- variable  model,  engineers  will  be  able  to  control  all  of  the  variables 
in  experimentation.  They  will  be  able  to  set  the  frequency  of  excitation  (ft))  by  adjusting 
the  environmental  conditions  affecting  the  piezoelectric  generator.  Additionally,  they  will 
be  able  to  set  the  natural  frequency  (ft)„)  and  modal  damping  ratio  (£)  of  the  generator  by 
adjusting  the  design  of  the  prototype. 

In  the  practical  application  of  a  piezoelectric  generator,  certain  conditions  may  not  be  as 
easy  to  control.  A  piezoelectric  generator  designed  to  operate  using  the  vibrations  of  an 
air  conditioning  vent  will  not  be  able  to  depend  on  a  precise  or  consistent  frequency  of 
excitation.  Engineers  must  design  the  generator  to  operate  within  a  range  of  frequency 
while  producing  the  necessary  power. 

Earlier  in  this  chapter,  we  determined  that  the  optimal  solution  to  the  three-variable  model 
occurred  when  ft)  =  ft),,.  Given  the  fact  the  environment  determines  the  frequency  of  exci¬ 
tation,  we  fixed  co  during  the  robustness  analysis  while  varying  co„  and  £. 

Single  Variable  Manipulation 

We  have  determined  the  optimal  solution  of  the  three-variable  model  depends  on  minimiz¬ 
ing  £  while  maximizing  ft)  and  setting  co  =  ft),,.  Due  to  this  knowledge  we  evaluate  the 
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robustness  of  the  optimal  solution  to  the  three-variable  model  in  terms  of  varying  £  and 
varying  the  difference  between  ft)  and  ft),,. 

First,  we  set  co  =  ft),,  =  3607T  and  vary  £  in  order  to  evaluate  the  power  lost  while  deviating 
from  the  optimal  value  of  £.  See  the  results  listed  in  Table  5.21. 


Table  5.21:  Single  variable  manipulation  of  £. 


c 

%  change  in  £ 

Power 

%  Loss 

0.005 

baseline 

7.23  x  1010 

baseline 

0.00505 

1% 

7.162  x  1010 

0.99% 

0.00525 

5% 

6.889  x  1010 

4.76% 

0.0055 

10% 

6.576  x  1010 

9.09% 

0.006 

20% 

6.028  x  1010 

16.67% 

0.01 

100% 

3.617  x  1010 

50% 

0.02 

300% 

1.808  x  1010 

74.75% 

A  1%,  5%,  or  10%  change  in  £  results  in  power  losses  of  less  than  10%.  A  change  as 
drastic  as  doubling  the  desired  value  of  £  results  in  only  a  50%  loss  in  power. 

Next,  we  set  £  =  0.005,  ft)  =  360k  and  vary  ft),,  to  evaluate  the  resulting  power  lost.  See 
the  results  in  Table  5.22. 
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Table  5.22:  Single  variab 

le  manipulation  < 

Of  (On . 

(On 

%  change  in  con 

Power 

%  Loss 

360k 

baseline 

7.23  x  1010 

baseline 

360.36 n 

0.1% 

6.948  x  1010 

3.94% 

359.64 k 

-0.1% 

6.962  x  1010 

3.75% 

360.9 K 

0.25% 

5.775  x  1010 

20.16% 

359. 1 7T 

-0.25% 

5.798  x  1010 

19.84% 

361.87T 

0.5% 

3.608  x  1010 

50.12% 

358.27 r 

-0.5% 

3.626  x  1010 

49.87% 

363.6 n 

1% 

I.444  x  1010 

80.04% 

356.47T 

-1% 

1.449  x  1010 

79.96% 

318k 

5% 

7.154  x  108 

99.01% 

342 k 

-5% 

7.160  x  108 

99.01% 

The  reader  will  notice  the  drastic  power  loss  when  the  separation  between  co  and  con  in¬ 
creases.  We  find  that  even  a  1%  discrepancy  between  co  and  con  results  in  a  power  loss 
greater  than  a  300%  change  in  £  from  its  optimal  value.  Whether  the  con  is  greater  than  or 
less  than  co  does  not  make  a  difference. 

Another  way  to  visualize  the  sensitivity  of  con  and  £  is  to  plot  the  resulting  power  in  terms 
of  optimal  power  as  each  variable  deviates  from  its  optimal  value.  See  Figure  5.2  below. 
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Figure  5.2:  Single  variable  manipulation  and  resulting  percent  optimal  power. 


As  seen  in  Figure  5.2,  the  reduction  of  power  when  £  increases  from  its  optimal  value  is 
linear,  with  a  relatively  small  slope.  On  the  other  hand,  when  ft)„  increases,  resulting  in 
an  increase  in  the  distance  between  co  and  con,  power  drastically  reduces  within  the  first 
percentages  of  change  in  to,,. 

Alternatively,  as  £  decreases  from  its  optimal  value,  power  increases  linearly  with  the  same 
small  slope.  The  reader  should  note,  reducing  £  from  its  optimal  value  sets  £’s  value  out¬ 
side  the  values  originally  studied  in  this  paper.  As  ft),,  decreases  from  its  optimal  value, 
resulting  in  an  increase  in  the  distance  between  co  and  ft)„,  power  drastically  reduces  sym¬ 
metrically  to  the  power  reduction  resulting  from  an  increase  in  ft),,  from  its  optimal  value. 

Upon  completion  of  the  single  variable  manipulation  robustness  analysis  of  the  three- 
variable  model,  we  determine  the  importance  of  minimizing  the  discrepancy  between  co 
and  ft),,.  £  should  be  set  as  low  as  possible,  but  the  precision  of  setting  £  should  be  a 
second  priority  to  precisely  matching  ft)  to  ft),,. 
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Multiple  Variable  Manipulation 

Single  variable  manipulation  provides  insight  into  how  the  departure  from  the  optimal  value 
of  each  variable  affects  power  output  when  the  other  variables  are  fixed  at  their  optimal  val¬ 
ues.  While  single  variable  manipulation  provides  valuable  analysis,  it  is  seldom  realistic  to 
assume  fixing  variable  values  in  practical  application.  A  more  realistic  analysis  comprises 
of  treating  each  variable  value  as  a  function  of  a  distribution.  The  distribution  has  an  upper 
and  lower  limit  to  match  the  possible  limits  guaranteed  by  a  manufacture,  for  example,  a 
capacitor  providing  10±1  ohms  of  resistance.  We  then  analyze  how  the  variable  distri¬ 
butions  affect  the  mean  and  median  of  power  output.  From  this  information  we  have  the 
ability  to  make  informed  recommendations  in  allowing  more  or  less  variability  in  variable 
distributions. 

We  assume  that  engineers  would  design  the  piezoelectric  generator  in  order  to  optimize 
power  output  with  a  specific  frequency  of  excitation  in  mind.  The  generator  itself  cannot 
manipulate  the  environment  and  accompanying  frequency  of  excitation.  To  simulate  this 
assumption  we  set  the  frequency  of  excitation,  ft),  to  the  optimal  value  of  3607T.  We  are 
more  interested  in  finding  how  the  variation  in  design  specifications  of  the  piezoelectric 
generator  affects  power  output. 
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In  order  to  guarantee  that  variable  values  stayed  within  desired  limits,  we  use  triangle 
distributions.  First,  using  SAS’s  commercial  statistical  discovery  software,  JMP,  we  create 
triangle  distributions  for  both  ft),,  and  £  with  modes  equal  to  their  respective  optimal  values 
(ft)„= 3607T,  £=0.005)  and  the  limits  at  ±10%  the  optimal  values.  Using  JMP,  we  create 
10,000  sample  points  from  the  distributions  and  fixed  value  of  ft)  and  computed  10,000 
power  outputs  using  the  three-variable  model.  The  distributions  of  ft),, ,  £  and  the  resulting 
power  output  are  illustrated  in  Figure  5.3. 
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Figure  5.3:  ft),,  ±  10%,  £  ±  10%,  and  resulting  power  distributions. 


Letting  ft),,  and  £  vary  10%  from  their  optimal  values  results  in  a  mean  power  output  of 
9.976  x  109.  This  mean  power  represents  an  86%  loss  from  the  optimal  value.  Furthermore, 
the  median  power  output  is  only  2.09  x  109,  a  97%  loss  from  the  optimal  power  value.  This 
means  50%  of  all  the  design  points  produce  a  loss  of  97%  or  more. 
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Next,  we  use  insights  derived  from  the  analysis  of  the  single  variable  manipulation  to  drive 
our  quest  for  more  power  output.  We  see  from  single  variable  manipulation  the  importance 
of  setting  co  —  (On  in  order  to  maximize  power.  In  the  following  distribution  we  allow  (On 
to  vary  only  5%  from  its  optimal  value,  while  allowing  £  to  vary  10%.  See  Figure  5.4  for 
the  0)n,  £,  and  resulting  power  distributions. 
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Figure  5.4:  (On±  5%,  £±10%,  and  resulting  power  distributions. 


Tightening  the  distribution  to  allow  for  only  a  5%  variation  of  con  while  keeping  the  previ¬ 
ously  used  £  distribution  results  in  a  mean  power  output  of  1.782  x  1010,  a  75%  loss  from 
the  optimal  power  value.  The  resulting  median  power  output  for  these  distributions  equals 
7.51  x  109,  a  90%  power  loss.  Even  if  the  results  are  not  overwhelming,  this  restriction  of 
variation  for  (On  leads  to  an  improvement  in  power. 
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Due  to  the  results  we  see  from  slightly  tightening  the  ft),,  distribution,  we  decide  to  inves¬ 
tigate  further  by  allowing  ft),,  to  vary  only  1%  from  its  optimal  value  while  allowing  £  to 
continue  to  vary  10%  from  optimal.  See  Figure  5.5  for  the  ft),, ,  £,  and  resulting  power 
distributions. 
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Figure  5.5:  co„±l%,  ^ dr  10%,  and  resulting  power  distributions. 


Further  tightening  the  distribution  to  allow  for  only  a  1%  variation  of  ft),,  while  keeping 
the  previously  used  £  distribution  results  in  a  mean  power  output  of  5.131  x  1010,  a  29% 
loss  from  the  optimal  power  value.  The  resulting  median  power  output  for  these  distribu¬ 
tions  equals  5.4  x  1010,  a  25%  power  loss.  This  is  a  drastic  improvement  in  power  output. 
Notice  the  power  distribution’s  shape  changes  as  well,  resulting  in  an  increased  likelihood 
of  generating  power  close  to  the  optimal  and  shifting  the  median  power  output  above  the 
mean. 
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Next,  we  try  to  further  improve  upon  our  findings  by  restricting  the  variation  of  £  from 
its  optimal  value.  We  tighten  the  distribution  of  £  to  only  5%  of  its  optimal  value  while 
keeping  ft),,’ s  distribution  to  allow  an  deviation  of  1%  of  its  optimal  value.  See  Figure  5.6 
for  the  (on,  £,  and  resulting  power  distributions. 
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Figure  5.6:  C0„=tl%,  £±5%,  and  resulting  power  distributions. 


Tightening  the  distribution  to  allow  for  only  a  5%  variation  of  £  while  keeping  the  pre¬ 
viously  used  (on  distribution  results  in  a  mean  power  output  of  5.128  x  1010,  a  29%  loss 
from  the  optimal  power  value.  The  resulting  median  power  output  for  these  distributions 
equals  the  median  of  the  previous  distributions  at  5.4  x  1010,  a  25%  power  loss.  We  find 
tightening  £’s  distribution  actually  leads  to  a  decrease  in  mean  power  output.  This  is  due  to 
the  fact  when  ft)  =  ft),,  as  £  approaches  zero  power  increases  infinitely.  That  is,  there  exists 
a  vertical  asymptote  when  £  =  0.  With  a  tighter  distribution,  the  value  of  £  does  not  have 
the  ability  to  take  on  smaller  values. 

Using  multiple  variable  manipulation  robust  analysis  we  find  when  attempting  to  optimize 
power  it  is  best  to  allow  as  little  variation  as  possible  from  natural  frequency’s  (ft),,)  optimal 
value.  This  will  minimize  the  discrepancy  between  ft)  and  ft),,  and  result  in  greater  power 
output.  The  restriction  of  the  variation  of  the  modal  damping  ration  (£)  from  its  optimal 
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value  is  not  as  important.  In  fact,  we  find  if  the  variation  of  £  would  allow  a  smaller  value 
than  the  desired  setting,  the  model  generates  even  greater  power. 


5.5.2  Six-Variable  Model 

In  the  six-variable  model,  engineers  are  able  to  control  all  of  the  variables  in  experimen¬ 
tation.  They  are  able  to  set  the  base  motion  driving  frequency  ((Dbase)  by  adjusting  the 
environmental  conditions  affecting  the  piezoelectric  generator.  Additionally,  they  are  able 
to  set  the  electromechanical  coupling  coefficient  (A),  load  resistance  (R),  net  clamped  ca¬ 
pacity  of  piezoelectric  material  (Co),  modal  short  circuit  net  frequency  ( wsc )  and  modal 
damping  ratio  (£)  by  adjusting  the  design  of  the  prototype. 

In  the  practical  application  of  a  piezoelectric  generator  certain  conditions  may  not  be  as 
easy  to  control.  A  piezoelectric  generator  designed  to  operate  using  the  vibrations  of  an  air 
conditioning  vent  will  not  be  able  to  depend  on  a  precise  or  consistent  base  motion  driving 
frequency.  Given  the  fact  the  environment  determines  the  base  motion  driving  frequency, 
we  fix  (ObaSe  during  the  robustness  analysis  while  varying  the  other  five  factors. 


Single  Variable  Manipulation 

Earlier  in  this  chapter,  we  determined  the  optimal  solutions  of  the  six-variable  model  in 
order  to  optimize  power.  To  remind  the  reader  of  these  values  we  include  Table  5.23. 


Table  5.23:  Optimal  values  for  each  variable  in  the  six-variable  model. 


Max  Power 

A 

R 

®base 

Co 

Wsc 

c 

1.8256  x  1010 

0.3376 

100 

3607 z 

1  x  10~6 

360 Ti 

0.005 

First,  we  fix  all  variables  at  their  optimal  values  and  vary  the  electromechanical  coupling 
coefficient  (A).  See  the  results  in  Table  5.24. 
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Table  5.24:  Single  variable  manipulation  of  A. 


A 

%  change  in  A 

Power 

%  Loss 

0.3376 

baseline 

1.8256  x  1010 

baseline 

0.30384 

-10% 

1.8057  x  1010 

1.09% 

0.32072 

-5% 

1.8209  x  1010 

0.26% 

0.3342 

-1% 

1.8254  x  1010 

0.01% 

0.3410 

1% 

1.8254  x  1010 

0.01% 

0.3545 

5% 

1.8211  x  1010 

0.25% 

0.3714 

10% 

1.8088  x  1010 

0.92% 

From  the  data  presented,  it  appears  varying  A  does  not  drastically  decrease  power  output. 
Varying  A  by  10%  from  its  optimal  value  only  results  in  around  a  1%  loss  in  power. 

Next,  we  fix  all  variables  at  their  optimal  values  and  vary  load  resistance  ( R ).  See  the  results 
in  Table  5.25. 


Table  5.25:  Single  variable  manipulation  off?. 


R 

%  change  in  R 

Power 

%  Loss 

100 

baseline 

1.8256  x  1010 

baseline 

90 

-10% 

1.8175  x  1010 

0.44% 

95 

-5% 

1.8228  x  1010 

0.15% 

99 

-1% 

1.8252  x  1010 

0.02% 

101 

1% 

1.8259  x  1010 

-0.02% 

105 

5% 

1.8262  x  1010 

-0.04% 

110 

10% 

1.8251  x  1010 

0.03% 

By  analyzing  the  data  above  we  see  that  varying  R  changes  power  by  even  less  than  the 
variation  of  A.  The  reader  will  notice  negative  loss  when  R  takes  on  the  value  of  101  and 
105  ohms.  This  represents  a  gain  of  power  at  those  values  of  R.  The  reader  will  also  notice 
increasing  R  past  those  values  to  110  ohms  results  in  a  power  loss  from  the  baseline.  It  is 
interesting  to  note,  that  had  we  initially  included  105  ohms  in  the  range  of  load  resistance 
to  be  studied  we  would  find  a  different  optimal  value  of  R. 
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Next,  we  fix  all  variables  at  their  optimal  values  and  vary  the  net  clamped  capacity  of 
piezoelectric  material  (Co).  See  the  results  in  Table  5.26. 


Table  5.26:  Single  variable  manipulation  of  Cq. 


C0 

%  change  in  Co 

Power 

%  Loss 

1  x  10~6 

baseline 

1.8256  x  1010 

baseline 

9  x  10~7 

-10% 

1.8223  x  1010 

0.18% 

9.5  x  10-7 

-5% 

1.8239  x  1010 

0.09% 

9.9  x  10~7 

-1% 

1.8253  x  1010 

0.02% 

1.01  x  10~6 

1% 

1.8259  x  1010 

-0.02% 

1.05  x  10“6 

5% 

1.8274  x  1010 

-0.1% 

1.1  x  10^6 

10% 

1.8292  x  1010 

-0.2% 

The  variation  of  Co  from  its  optimal  value  results  in  small  variation  of  power  output.  Also, 
as  Co  increases  so  does  the  power.  From  the  results  above,  we  assume  if  we  had  included 
a  larger  value  as  the  maximum  limit  of  Co  to  study  we  would  derive  a  larger  value  of  Co  as 
its  optimal  value. 

Next  we  fix  all  variables  at  their  optimal  values  and  vary  the  modal  short  circuit  net  fre¬ 
quency  (vf5c).  See  the  results  in  Table  5.27. 


Table  5.27:  Single  variable  manipulation  of  wsc- 


Wsc 

%  change  in  wsc 

Power 

%  Loss 

36071 

baseline 

1.8256  x  1010 

baseline 

324k 

-10% 

2.0229  x 108 

98.89% 

342 n 

-5% 

7.5408  x  108 

95.87% 

336.4k 

-1% 

9.7642  x  109 

46.51% 

363.6 K 

1% 

8.5501  x  109 

53.17% 

3787T 

5% 

6.5336  x  108 

96.42% 

396 n 

10% 

1.6194  x  108 

99.11% 

Variation  in  wsc  results  in  large  power  loss.  As  small  as  a  1%  deviation  in  modal  short 
circuit  frequency  results  in  close  to  a  50%  loss  in  power.  A  10%  deviation  of  wsc  results  in 
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a  99%  power  loss.  The  six-variable  model  is  extremely  sensitive  to  variations  of  wsc  from 
its  optimal  value. 

As  in  the  three-variable  model,  the  difference  between  two  frequencies  in  the  six-variable 
model  greatly  impact  power  generation.  In  this  robustness  study  we  fix  (O^ase  at  its  optimal 
value  and  vary  wsc ■  Due  to  the  fixed  COi)ase,  the  relationship  between  the  two  frequencies 
and  power  generation  is  not  obvious  to  the  reader.  However,  after  analyzing  Equation 
(4.1),  the  reader  will  notice  when  (di3ase  =  wsc  the  denominator  decreases  by  eliminating 
(  —  U>base2  +  W’sc2 )  ( 1  A  Qf  0)/7ase2 R2) .  The  decrease  in  the  denominator  of  Equation  (4.1) 
leads  to  greater  power  generation. 

Finally,  we  fix  all  variables  at  their  optimal  values  and  vary  the  modal  damping  ratio  (£). 
See  the  results  in  Table  5.28. 


Table  5.28:  Single  variable  manipulation  of  £. 


c 

%  change  in  £ 

Power 

%  Foss 

0.005 

baseline 

1.8256  x  1010 

baseline 

0.0045 

-10% 

2.0227  x  1010 

-10.79% 

0.00475 

-5% 

1.9203  x  1010 

-5.19% 

0.00495 

-1% 

1.8440  x  1010 

-1.01% 

0.00505 

1% 

1.8075  x  1010 

0.99% 

0.00525 

5% 

1.7377  x  1010 

4.82% 

0.0055 

10% 

1.6560  x  1010 

9.29% 

The  variation  of  £  from  its  optimal  value  results  in  small  variation  of  power  output.  Also, 
as  £  decreases  the  power  increases.  From  the  results  above  we  could  assume  that  if  we 
had  included  a  smaller  value  as  the  minimum  limit  of  £  to  study  we  would  derive  a  smaller 
value  of  £  as  its  optimal  value,  as  in  Figure  5.7. 
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Another  way  to  visualize  the  sensitivity  of  each  variable  is  to  plot  the  resulting  power  in 
terms  of  optimal  power  as  each  variable  deviates  from  its  optimal  value. 


Figure  5.7:  Single  variable  manipulation  and  resulting  percent  optimal  power. 


Figure  5.7  further  illustrates  the  need  to  set  wsc  to  its  precise  optimal  value.  The  effects  of 
all  other  variables  are  drowned  out  by  wsc’s  dominating  negative  effect.  Individually  A,  R, 
and  Co  change  power  output  insignificantly  even  under  ±10%  deviations  from  their  optimal 
values.  £  demonstrates  its  ability  to  linearly  effect  power  both  negatively  and  positively. 
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In  order  to  see  the  effects  of  the  less  sensitive  variables,  see  the  zoomed  in  plot  in  Figure 
5.8. 


Figure  5.8:  Zoomed-in  single  variable  manipulation  and  resulting  percent  optimal  power. 


From  Figure  5.8  one  will  notice,  after  wsc,  the  deviations  of  £  have  the  largest  impact  on 
power  output.  Also  notice  a  change  of  -10%  in  £  increases  power  slightly  more  than  a 
change  of  10%  in  £  decreases  power. 

The  zoomed-in  plot  also  illustrates  the  curvature  of  A’s  power  curve.  The  optimal  value  of 
A  would  not  have  changed  even  if  we  selected  a  wider  range  of  A  values  to  study. 

Upon  completion  of  the  single  variable  manipulation  of  the  six- variable  model  robustness 
analysis,  we  find  modal  short  circuit  net  frequency  (wsc)  proves  by  far  the  most  sensitive 
variable.  A  deviation  of  ±10%  wsc  from  its  optimal  value  results  in  around  a  99%  power 
loss.  On  the  opposite  end  of  the  spectrum  a  deviation  of  ±10%  A,  R,  or  Co  from  their 
optimal  values  results  in  a  power  loss  of  around  1%  or  less.  Interesting  to  note,  slightly 
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increasing  R  or  Co  or  decreasing  £  results  in  a  power  gain. 

Multiple  Variable  Manipulation 

As  previously  discussed  in  the  three-variable  model  robustness  analysis,  the  single  variable 
manipulation  of  the  six-variable  model  provides  insight  into  how  the  departure  from  the 
optimal  value  of  each  variable  affects  power  output  when  the  other  variables  are  fixed  at 
their  optimal  values.  Like  the  three-variable  multiple  variable  robustness  analysis  we  use 
variable  distributions  to  attempt  to  identify  which  variables  require  more  precision  and  thus 
less  variation  from  their  optimal  value  in  order  to  maximize  power. 

We  assume  that  engineers  design  the  piezoelectric  generator  in  order  to  optimize  power  out¬ 
put  with  a  specific  frequency  of  excitation  in  mind.  The  generator  itself  cannot  manipulate 
the  environment  and  accompanying  frequency  of  excitation.  To  simulate  this  assumption 
we  set  the  base  motion  driving  frequency,  (Qbase,  to  the  optimal  value  of  3607T.  We  are 
more  interested  in  finding  how  the  variation  in  design  specifications  of  the  piezoelectric 
generator  affect  power  output. 
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In  order  to  guarantee  that  variable  values  stay  within  desired  limits  we  use  triangle  distribu¬ 
tions.  First,  using  SAS’s  commercial  statistical  discovery  software,  JMP,  we  create  triangle 
distributions  for  the  five  variables  in  the  six-variable  model  with  modes  equal  to  their  re¬ 
spective  optimal  values  (A  —  0.3376,  R  —  100,  Co  =  1  x  10  6,  wsc  —  360 n,  and  £=0.005) 
and  the  limits  at  ±  10%  the  optimal  values.  Using  JMP,  we  create  10,000  sample  points 
from  the  distributions  and  fixed  value  of  (0/yase,  and  compute  10,000  power  outputs  using 
the  six-variable  model.  See  the  distributions  of  the  variables  and  the  resulting  power  output 
in  Figure  5.9. 
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Figure  5.9:  A±  10%,  R±  10%,  Co±  10%,  wsc±  10%,  £ ±  10%,  and  resulting  power  distributions. 


Letting  the  five  variables  vary  10%  from  their  optimal  values  results  in  a  mean  power 
output  of  4.4875  x  109.  This  mean  power  represents  a  75%  loss  from  the  optimal  value. 
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Furthermore,  the  median  power  output  is  only  1.92  x  109,  an  89%  loss  from  the  optimal 
power  value.  This  means  50%  of  all  the  design  points  produce  a  loss  of  89%  or  more  of 
the  optimal  value. 


Next,  we  use  insights  derived  from  the  analysis  of  the  single  variable  manipulation  to  seek 
more  power  output.  From  single  variable  manipulation  we  determine  wsc  to  be  the  most 
sensitive  variable  in  the  six-variable  power  generation  model.  In  the  following  distribution 
we  allow  wsc  to  vary  only  5%  from  its  optimal  value,  while  allowing  the  other  four  vari¬ 
ables  a  10%  variation.  See  the  distributions  of  the  variables  and  the  resulting  power  output 
in  Figure  5.10. 
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Figure  5.10:  A±  10%,  R±  10%,  Cq±  10%,  wsc± 5%,  £ ±  10%,  and  resulting  power  distributions. 


Allowing  wsc  to  vary  only  5%  from  its  optimal  value  while  allowing  the  other  four  vari- 
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ables  a  10%  variation  results  in  a  mean  power  of  7.6514  x  109,  representing  a  58%  power 
loss.  This  is  a  17%  improvement  in  mean  power  loss  from  allowing  wsc  10%  variability. 
Additionally,  the  median  power  output  is  5.8  x  109,  a  68%  power  loss.  The  median  power 
loss  is  a  21%  improvement  over  allowing  wsc  10%  variability. 


Figure  5.10  illustrates  the  benefit  of  constraining  the  variation  of  the  modal  short  circuit 
net  frequency.  Next,  we  further  tighten  the  variation  of  wsc  to  1%  while  maintaining  the 
other  four  variables  at  10%  variation.  See  the  distributions  of  the  variables  and  the  resulting 
power  output  in  Figure  5.11. 
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Figure  5.11:  A±  10%,  R±  10%,  Cq±  10%,  wsc±  1%,  £ ±  10%,  and  resulting  power  distributions. 


Further  tightening  the  variation  of  the  modal  short  circuit  net  frequency  results  in  an  even 
greater  mean  power  output  of  1.599  x  1010,  representing  only  a  12%  power  loss  from  the 
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optimal  power.  Furthermore,  the  reader  will  notice  the  power  distribution  changing  drasti¬ 
cally  as  well,  resulting  in  a  median  power  greater  than  the  mean  power.  The  median  power 
output  is  1.7  x  1010,  a  power  loss  of  only  7%  from  the  optimal  power. 


We  learn  the  importance  of  minimizing  the  variation  of  the  modal  short  circuit  frequency 
from  its  optimal  value.  If  we  cannot  guarantee  a  1%  variation  constraint  on  wsc ,  is  it 
possible  to  make  up  power  output  by  tightening  the  variation  of  the  other  four  variables? 
We  attempt  to  answer  the  above  question  by  allowing  all  five  variables  a  5%  variation  to 
compare  power  output  to  previous  results.  See  the  distributions  of  the  variables  and  the 
resulting  power  output  in  Figure  5.12. 
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Figure  5.12:  A  ±5%,  R±  5%,  Co  ±5%,  wsc±5%,  £±5%,  and  resulting  power  distributions. 


From  Figure  5.12  we  see  a  power  distribution  similar  to  that  of  Figure  5.10.  In  fact,  restrict- 
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ing  all  variables  to  a  5%  variation  results  in  a  mean  power  of  7.656  x  109,  representing  a 
58%  power  loss  from  the  optimal  power.  This  mean  power  is  extremely  similar  to  the  re¬ 
sults  displayed  in  Figure  5.10.  Additionally  similar  to  the  results  from  Figure  5.10,  the 
median  power  output  is  5.75  x  109,  representing  a  69%  power  loss.  We  find  tightening  the 
variation  of  the  variables  other  than  the  modal  short  circuit  net  frequency  results  in  negligi¬ 
ble  power  improvement  over  when  we  restrict  wsc  to  5%  variation  and  allow  10%  variation 
for  the  other  four  variables. 


Our  next  question  asks  whether  it  is  possible  to  generate  a  best  case  scenario  by  restricting 
all  variables  to  1%  variation.  Do  we  see  significant  improvement  over  when  we  restrict  wsc 
to  1%  variation  while  allowing  the  four  other  variables  10%  variation?  See  the  distributions 
of  the  variables  and  the  resulting  power  output  in  Figure  5.13. 
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Figure  5.13:  A  ±  1%,  R±  1%,  Cq ±  1%,  wsc  ±  1%,  C  3=  1%,  and  resulting  power  distributions. 
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Comparing  the  power  distribution  in  Figure  5.13  with  the  power  distribution  in  Figure  5.11, 
we  notice  the  likelihood  of  producing  power  close  to  the  optimal  increases.  The  power 
distribution  in  Figure  5.11,  however,  displays  the  possibility  of  producing  power  greater 
than  the  optimal.  As  we  discussed  in  single  variable  manipulation,  this  occurs  when  £  is 
smaller  than  its  optimal  value  or  Co  or  R  is  greater  than  its  optimal  value. 

Restraining  all  variables  to  a  1%  variation  results  in  a  mean  power  of  1.601  x  1010,  repre¬ 
senting  a  12%  power  loss,  and  a  median  power  of  1.7  x  1010,  representing  a  7%  power  loss. 
The  restriction  of  all  variables  to  a  1%  variation  results  in  a  negligible  benefit  in  terms  of 
mean  and  median  power  compared  to  restraining  modal  short  circuit  net  frequency  to  1% 
variation  while  allowing  the  other  four  variables  10%  variation. 

We  see  when  we  tighten  the  variation  of  wsc  to  1%  and  allow  10%  variation  of  the  other 
variables  there  is  a  possibility  of  generating  power  greater  than  the  optimal.  During  single 
variable  manipulation  we  determined  decreasing  £  from  its  optimal  value  has  the  greatest 
influence  in  increasing  power.  We  also  note  decreasing  £  10%  increases  power  slightly 
more  than  increasing  £  10%  decreases  power.  Are  we  able  to  produce  a  higher  mean 
and  median  power  if  we  allow  £  to  vary  10%  while  restricting  all  other  variables  to  1% 
variation?  See  the  distributions  of  the  variables  and  the  resulting  power  output  in  Figure 
5.14. 
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Figure  5.14:  A±  1%,  R±  1%,  Cq±  1%,  1%.  C  ±  10%,  and  resulting  power  distributions. 


The  power  distribution  in  Figure  5.14  looks  extremely  similar  to  the  power  distribution  in 
Figure  5.11  when  we  restrict  wsc  to  1%  variation  while  allowing  10%  variation  for  all  other 
variables.  For  both  power  distributions  there  are  possibilities  of  producing  power  greater 
than  the  optimal  and  a  high  likelihood  of  producing  power  approximate  to  the  optimal. 
Additionally,  the  power  distributions  have  extremely  similar  means  and  medians.  The  mean 
and  median  power  produced  for  the  distribution  above  are  1.604  x  1010  and  1.7  x  1010, 
respectively.  This  accounts  for  a  mean  power  loss  of  12%  and  a  median  power  loss  of  7%. 

Although  not  an  overwhelming  front  runner,  we  find  that  the  best  mean  power  results  from 
restricting  all  variables  to  1%  variation  while  allowing  £  10%  variation.  This  produces  a 
difference  of  only  3  x  107  greater  power  than  restricting  all  variables  to  1%  variation. 
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The  biggest  takeaway  from  the  six-variable  robustness  analysis  is  the  importance  in  mini¬ 
mizing  the  variation  of  the  modal  short  circuit  net  frequency  (wsc)  from  its  optimal  value. 
The  restriction  of  variation  for  all  other  variables  should  be  a  distant  second  priority.  In 
fact  allowing  the  modal  damping  ratio  (£)  to  vary  from  its  optimal  improves  mean  power 
production  slightly. 
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CHAPTER  6: 
Conclusion 


6.1  Summary  of  Results 

The  military  is  looking  for  alternative  energy  sources  to  provide  deployed  bases  with  elec¬ 
tricity.  Diesel  fueled  generators  and  their  large  logistics  tails  have  proven  expensive  and  a 
liability  to  force  protection.  The  alternative  energy  focus  for  now  is  solar,  but  solar  arrays 
require  infrastructure  and  space.  Piezoelectric  generators,  combined  with  other  renewable 
energy  means,  could  provide  a  more  efficient  means  of  electricity  production  to  fossil  fuel 
powered  generators.  Additionally,  new  methods  of  energy  storage  could  eventually  facili¬ 
tate  alternative  energy  use  in  military  applications  [33]. 

This  paper  attempted  to  find  the  optimal  designs  involving  material  parameters  in  the  piezo¬ 
electric  generator  as  well  as  the  generator’s  environment  in  order  to  maximize  electric  out¬ 
put.  Our  hope  is  that  the  results  of  this  paper  will  serve  as  a  step  towards  the  practical 
application  of  efficient  piezoelectric  generators  and  contribute  to  preserving  the  combat  ca¬ 
pability  of  the  US  Military.  Through  the  study  of  two  mathematical  models  we  found  the 
optimal  values  for  variables  given  constrained  ranges.  We  presented  the  optimal  values  for 
the  three-variable  model  in  Table  5.1  and  the  optimal  values  for  the  six-variable  model  in 
Table  5.12. 

6.1.1  Three-Variable  Model  Analysis 

In  the  detailed  analysis  of  the  three-variable  model,  we  study  the  power  as  both  a  function  of 
normalized  frequency  and  a  function  of  modal  damping  ratio  and  frequency.  When  we  ana¬ 
lyze  power  as  a  function  of  normalized  frequency,  we  find  that  generating  maximum  power 
the  optimal  normalized  frequency  (x  =  depends  on  the  value  of  the  modal  damping 
ratio  (£).  As  £  approaches  zero  the  optimal  normalized  frequency  approaches  one.  When 
we  analyze  power  as  a  function  of  modal  damping  ratio  and  normalized  frequency,  we  are 
able  to  visualize  the  potential  power  generation  capability  by  minimizing  £  and  setting  the 
normalized  frequency  to  one.  Finally,  we  use  Roundy  and  Wright’s  results  [15]  to  support 
our  simple  three- variable  model. 
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6.1.2  Six-Variable  Model  Analysis 

Through  the  analysis  of  the  six-variable  model  optimization  results  we  find  that  power  in¬ 
creases  as  the  base  motion  driving  frequency  (C0/,aye),  net  clamped  capacity  of  piezoelectric 
material  (Co),  and  modal  short  circuit  net  frequency  (wsc)  increases  and  modal  damping 
ratio  (£)  decreases.  The  optimal  values  for  the  electromechanical  coupling  coefficient  (A) 
and  load  resistance  ( R )  vary  depending  on  sampling  method.  There  appears  to  be  a  corre¬ 
lation  in  power  output  and  the  difference  of  CO^ase  and  wsc- 

6.1.3  Sampling  Methods 

In  our  search  for  optimal  values  to  maximize  power  production  we  use  two  sampling  meth¬ 
ods,  mr  factorial  and  nearly  orthogonal  and  balance  Latin  hypercubes  (NOBLH).  Although 
the  NOBLH  design  provides  nearly  orthogonal  samplings  with  fair  space-filling  properties, 
its  failure  to  effectively  sample  the  extremes  of  the  design  space  leads  to  less  than  optimal 
results.  In  both  of  our  models,  the  optimal  variable  values  exists  at  the  extremes  of  the 
variable  ranges  in  all  but  one  variable  (A).  mK  variable  sampling  chooses  design  points 
incorporating  the  optimal  values,  while  NOBLH  does  not. 

Had  the  models  required  more  computational  rigor  (more  factors  or  computations)  the 
NOBLH  samplings’  efficiency  might  have  provided  insights  into  optimization  where  mr 
factorial  design  would  have  been  too  costly.  However,  with  both  our  three-  and  six-variable 
models  we  do  not  have  difficulty  evaluating  mk  factorial  samplings  even  with  double  digit 
levels  (m). 

6.1.4  Iterative  Method 

The  NOBLH  iterative  methods  we  use  for  the  three-  and  six-variable  models  are  ad  hoc 
methods  that  arose  from  trying  to  glean  more  insight  about  our  models’  behavior.  Had 
we  been  interested  solely  in  identifying  an  optimal  solution,  then  an  adaptive  sequential 
procedure,  such  as  response  surface  methodology  [34],  STRONG  [35],  or  other  simulation 
optimization  approaches  would  require  much  less  sampling. 

6.1.5  Robustness  Analysis 

After  determining  the  optimal  variable  values,  our  analysis  turns  to  evaluating  the  robust¬ 
ness  of  the  variables  by  evaluating  power  generation  while  incorporating  deviations  from 
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the  optimal  variable  values.  In  the  three- variable  model,  we  find  the  separation  between  the 
frequency  of  excitation  of  the  base  (ft))  and  natural  frequency  of  the  generator  (ft),,)  to  be 
the  most  sensitive  variable.  While  the  designers  of  the  piezoelectric  generator  should  seek 
to  minimize  the  generator’s  modal  damping  ratio  (£),  the  designer’s  priority  should  be  to 
precisely  match  the  generator’s  natural  frequency  with  the  frequency  of  excitation  driving 
the  generator. 

The  precise  matching  of  the  base  motion  driving  frequency  {(Obase)  to  the  modal  short  circuit 
net  frequency  (wsc)  proves  to  be  the  most  important  factor  in  minimizing  power  loss  during 
the  robustness  analysis  of  the  six-variable  model.  Setting  the  electromagnetic  coupling 
coefficient  (A),  load  resistance  (R),  and  net  clamped  capacity  of  piezoelectric  material  (Co) 
to  their  respective  optimal  value  is  a  much  lower  priority  for  the  piezoelectric  generator 
designer.  In  fact,  our  robustness  study  concludes  if  the  likelihood  of  negative  variation  is 
greater  than  or  equal  to  positive  variation  for  the  design’s  modal  damping  ratio,  it  is  best 
to  let  £  vary  from  its  optimal  value  while  restraining  the  variation  of  other  variables  from 
their  optimal  values  in  order  to  produce  the  greatest  mean  power  output. 

6.2  Further  Studies 

Through  the  analysis  of  two  mathematical  models,  we  determine  optimal  values  for  each 
variable,  given  a  practical  variable  range  in  order  to  maximize  power  generation  of  a  piezo¬ 
electric  generator.  Additional  studies  are  needed  to  propel  our  findings  into  practical  appli¬ 
cation  of  piezoelectric  materials  for  energy  harvesting. 

6.2.1  PDE  Solution 

To  derive  the  nine  variable  model  in  Chapter  3,  we  use  a  one-mode  expansion  method. 
An  alternative  methods  to  derive  models  would  be  to  use  a  partial  differential  equation 
expansion  [36].  This  mathematical  derivation  would  provide  an  additional  piezoelectrical 
generator  model  to  conduct  further  optimization  studies. 

6.2.2  Variable  Range  Adjustment 

We  derive  the  variable  ranges  through  research  of  materials  specifications  and  previous 
experimental  measurements.  Personnel  more  experienced  in  vibration,  piezoelectric  mate¬ 
rials,  and  electrical  engineering  will  have  the  ability  to  further  adjust  the  ranges  for  each 
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variable.  Using  our  methods  of  analysis  and  power  generation  models,  it  is  possible  to  find 
updated  optimal  values  for  each  variable. 

6.2.3  Simulation  Optimization 

In  future  studies,  an  automated  simulation  optimization  technique  would  be  a  better  ap¬ 
proach  for  finding  optimal  nominal  settings  for  the  variables  while  keeping  sampling  re¬ 
quirements  to  a  minimum. 

6.2.4  Robustness  Analysis  Expansion 

Our  efforts  study  the  robustness  of  the  two  power  generation  models  based  on  both  single 
and  multiple  variable  manipulation.  In  both  the  three-  and  six-variable  power  manipulation, 
we  assume  the  environment  to  be  static.  In  other  words,  we  assume  the  frequency  of 
excitation  (to)  and  the  base  motion  driving  frequency  to  be  constant  for  the  three-  and 
six-variable  models,  respectively.  This  enables  us  to  focus  our  robustness  study  on  fewer 
variables. 

In  practical  applications,  the  frequency  of  the  driving  vibration  may  not  prove  constant. 
The  vibration  of  an  air  conditioning  duct  may  change  based  on  temperature,  humidity,  fan 
speed,  or  a  myriad  other  factors.  Further  robustness  analysis  focused  on  the  base  motion 
driving  frequency  may  lead  to  a  better  piezoelectric  generator  power  output  in  varying 
environmental  conditions. 

6.2.5  Model  Adjustment 

The  two  mathematical  models  we  use  do  not  serve  as  the  only  models  to  simulate  piezoelec¬ 
tric  power  production.  Further  studies  could  include  additional  models  or  adjustments  to 
the  two  models  used  in  this  study.  The  optimization  and  robustness  analysis  methodologies 
in  this  study  could  serve  as  a  framework  in  the  study  of  new  models. 

6.2.6  Practical  Experimentation 

We  gear  the  efforts  of  this  paper  in  order  to  find  insights  into  maximizing  piezoelectric 
generator  power  output.  Using  our  findings,  engineers  can  attempt  to  match  our  modeled 
results  using  a  prototype  piezoelectric  generator  designed  with  our  optimal  variable  val¬ 
ues.  Feedback  from  the  results  of  the  experimentation  could  serve  to  adjust  our  models  or 
variable  limits. 
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APPENDIX:  MATLAB  Codes 


A.l  Three- Variable  Orthogonality  Computation  Code 

%3  Variable  Orthogonality  Computation 
%Created  April  2015  by  Russell  Nelson 

clear  all 

%Import  design  points 
filename  =  ' Factorial3var . csv ' ; 

%Rename  file 
NOLH=csvread (filename) ; 

%Compute  orthogonality  for  each  variable 

avgomega  =  mean (NOLH ( : , 1 ) ) ; 

avgomega_n  =  mean (NOLH ( : , 2 ) ) ; 

avgzeta  =  mean (NOLH ( :  ,  3 ) )  ; 

numeratorl  =  0; 

numerator2  =  0; 

numerator3  =  0; 

denomll  =  0; 

denoml2  =  0; 

denom21  =  0; 

denom22  =  0; 

denom31  =  0; 

denom32  =  0; 

for  i  =  1:1000 

omegai  =  NOLH(i,l); 
omega_ni  =  NOLH  (i, 2); 

numeratorl  =  numeratorl  +  (omegai-avgomega) * (omega_ni-avgomega_n) ; 
denomll  =  denomll  +  (omegai-avgomega) ^2 ; 
denoml2  =  denoml2  +  (omega_ni-avgomega_n) A2; 

end 

%Outputs  correlation  for  omega  to  omega_n 
Ortho_omega_omega_n=numeratorl/ sqrt (denoml l*denoml2 ) 
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for  i  =  1:1000 

omega!  =  NOLH(i,l); 
zetai  =  NOLH  (i, 3) ; 

numerator!  =  numerator!  +  (omegai-avgomega) * (zetai-avgzeta) ; 
denom21  =  denom21  +  (omegai-avgomega) ^2 ; 
denom22  =  denom22  +  ( zetai-avgzeta) ^2 ; 

end 

%Outputs  correlation  for  omega  to  zeta 
Ortho_omega_zeta  =  numerator2 /sqrt (denom2 1 *denom22 ) 

for  i  =  1:1000 

omega_ni  =  NOLH (i, 2); 
zetai  =  NOLH  (i, 3) ; 

numerator3  =  numerator3  +  (omega_ni-avgomega_n) *( zetai-avgzeta) 
denom31  =  denom31  +  (omega_ni-avgomega_n) A2; 
denom32  =  denom32  +  ( zetai-avgzeta) ^2 ; 

end 

%Outputs  correlation  for  omega_n  to  zeta 
Ortho_omega_n_zeta  =  numerator3/sqrt (denom31*denom32) 


A.2  Six- Variable  Orthogonality  Computation  Code 

%6  Variable  Orthogonality  Computation 
%Created  April  2015  by  Russell  Nelson 

clear  all 

%Import  design  points 
filename  =  ' Factorial6var . csv ' ; 

%Rename  file 
NOLH=csvread (filename) ; 

%Compute  orthogonality  for  each  variable 
num=0 ; 
denl=0 ; 
den2=0 ; 


84 


for  d  =  2 : 6 

for  i  =1:15625 

Ai  =  NOLH  ( i , 1)  ; 

Bi  =  NOLH (i, d) ; 

avA  =  mean (NOLH ( : , 1 ) ) ; 

avB  =  mean (NOLH ( : , d) ) ; 

nun  =  num  +  (Ai-avA) * (Bi-avB) 

deni  =  deni  + (Ai-avA) A2; 

den2  =  den2  + (Bi-avB) A2; 

end 

OrthoA(d)=  num/sqrt (denl*den2 ) ; 

end 

num=0  ; 
denl=0 ; 
den2=0 ; 
for  d  =  3:6 

for  i  =1:15625 

Ai  =  NOLH  (i, 2)  ; 

Bi  =  NOLH (i, d) ; 

avA  =  mean (NOLH  (  : , 2) ) ; 

avB  =  mean (NOLH ( : , d) ) ; 

num  =  num  +  (Ai-avA) * (Bi-avB) 

deni  =  deni  + (Ai-avA) A2; 

den2  =  den2  + (Bi-avB) A2; 

end 

OrthoR(d)=  num/sqrt (denl*den2 ) ; 

end 

num=0  ; 
denl=0 ; 
den2=0 ; 
for  d  =  4:6 

for  i  =1:15625 

Ai  =  NOLH (i, 3) ; 

Bi  =  NOLH (i, d) ; 

avA  =  mean (NOLH ( : , 3) ) ; 

avB  =  mean (NOLH ( : , d) ) ; 

num  =  num  +  (Ai-avA) * (Bi-avB) 
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deni  =  deni  + (Ai-avA) A2; 
den2  =  den2  + (Bi-avB) A2; 

end 

OrthoO(d)=  num/sqrt (denl*den2 )  ; 

end 

num=0  ; 
denl=0 ; 
den2=0 ; 
for  d  =  5:6 

for  i  =1:15625 

Ai  =  NOLH  ( i , 4)  ; 

Bi  =  NOLH (i, d) ; 

avA  =  mean (NOLH  (  : , 4 ) ) ; 

avB  =  mean (NOLH ( : , d) ) ; 

nun  =  num  +  (Ai-avA) * (Bi-avB) ; 

deni  =  deni  + (Ai-avA) A2; 

den2  =  den2  + (Bi-avB) A2; 

end 

OrthoC(d)=  num/sqrt (denl*den2 )  ; 

end 

num=0 ; 
denl=0 ; 
den2=0 ; 
for  d  =  6:6 

for  i  =1:15625 

Ai  =  NOLH (i, 5) ; 

Bi  =  NOLH (i, d) ; 

avA  =  mean (NOLH ( : , 5) ) ; 

avB  =  mean (NOLH ( : , d) ) ; 

num  =  num  +  (Ai-avA) * (Bi-avB)  ; 

deni  =  deni  + (Ai-avA) A2; 

den2  =  den2  + (Bi-avB) A2; 

end 

OrthoW(d)=  num/sqrt (denl*den2 )  ; 

end 
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A.3  Euclidean  Maximin  Distance  Code 


%Euclidean  Maximin  Distance  Calculator 
%Created  April  2015  by  Russell  Nelson 
clear  all 

%Import  factorial  design  points  to  use 
%as  reference 

file  =  ' Factorial3var . csv ' ; 
base  =  csvread (file) ; 

%Import  Design  Points 
filename  =  ' Factorial3var . csv ' ; 

%Rename  file 
A=csvread (filename)  ; 

%Establish  number  of  variables  (n) 

%and  number  of  design  points  (k) 

[n,k]  =  size  (NOLH) ; 

d=0 ; 

i=l; 

%Normalizing  vectors 
for  p  =  1 : k 

avg  =  (max (base ( : , p) ) -min (base ( : , p) ) ) /2+min (base ( : , p) ) 
for  q  =  1 : n 

A ( q, p )  =  (A (q, p) -avg) /avg; 

end 

end 

%Computing  all  distances 
for  s  =  1 : n 

for  t  =  2 : n 

if  s  ==  t 
else 

for  j  =  1 : k 

d  =  d  +  abs (A ( s , j ) -A (t ,  j ) ) A2  ; 
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end 

d  =  sqrt (d) ; 
distance  (i)  =  d; 
d=0  ; 

i  =  i+1; 

end 

end 

end 

%Finding  the  Euclidean  maximin  distance 
min (distance ) 


A.4  Three-Variable  Model  Optimization  Code 
for  Factorial  Sampling 

%Optimal  Power  for  Three  Variable  Model  (Factorial  Sampling  Method) 
%Created  January  2015  by  Russell  Nelson. 

%The  following  script  evaluates  the  average  power  dissipated  by  the  load 
%resistor  of  a  piezoelectric  generator.  The  equation,  consisting  of  3 
%independent  variables,  was  derived  by  Professor  Hong  Zhou, 
clear  all 

a=l; 
b=l; 
c=l  ; 

%Set  variable  limits  (range) 
omega_min_lim  =  120*pi; 
omega_max_lim  =  360*pi; 
omega_n_min_lim  =  120*pi; 
omega_n_max_lim  =  360*pi; 
zeta_min_lim  =  0.005; 
zeta_max_lim  =  0.02; 

%Set  number  of  levels  per  variables 
levels  =  100; 
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%Computes  values  of  power  for  all  combinations  of 
%variables  at  all  levels 

for  omega  =  linspace (omega_min_lim, omega_max_lim, levels) 
b=l  ; 

for  omega_n  =  linspace (omega_n_min_lim, omega_n_max_lim, levels) 
c=l ; 

for  zeta  =  linspace ( zeta_min_lim, zeta_max_lim, levels ) 
P(a,b,c)  =  (zeta* ( (omega/omega_n) A3) * (omegaA3) ) / .  .  . 

(  ( 1-  (  (omegaA2 )  /  (omega_nA2 )  )  )  A2+  .  .  . 

(2*zeta* (omega/omega_n) ) A2)  ; 
c=c+l ; 

end 
b=b+l ; 

end 
a=a+l ; 

end 

%Vectorizes  P  matrix 
Pvec=P ( : )  ; 

%Find  the  maximum  power  in  Pvec 
[v, I ] =max (Pvec) 

%Determines  the  location  of  each  variable  level  for  maximum  power 
[d,  e, f ] =ind2sub (size(P),I); 

%Determines  the  value  for  each  variable  at  maximum  power 
omega  =  linspace (omega_min_lim, omega_max_lim, levels ) ; 
omegaopt=omega (d) 

omega_n  =  linspace (omega_n_min_lim, ome ga_n_ma x_l im , levels) ; 
omega_nopt=omega_n (e) 

zeta  =  linspace (zeta_min_lim, zeta_max_lim, levels) ; 
zetaopt=zeta (f ) 

A.5  Three-Variable  Model  Optimization  Code 
for  NOBLH  Sampling 
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%3  Variable  Model  Optimization  Code  for  NOBLH  Sampling 
%Created  January  2015  by  Russell  Nelson 

%The  following  script  evaluates  the  average  power  dissipated  by  the  load 
%resistor  of  a  piezoelectric  generator.  The  equation,  consisting  of  3 
%independent  variables,  was  derived  by  Professor  Hong  Zhou. 

clear  all 

%Import  design  points 
filename  =  ' NOL3var . csv ' ; 

%Rename  file 
NOLH=csvread (filename)  ; 

%Compute  power  for  all  design  points 
for  i  =1:512 

omega=NOLH  ( i , 1 ) ; 
omega_n=NOLH  (i, 2) ; 
zeta=NOLH  (i, 3) ; 

Power (l,i)  =  (zeta* ( (omega/omega_n) A3) * (omegaA3) ) / .  .  . 

( ( 1- ( ( omega A 2 ) / (omega_nA2 ) ) ) A2  + (2*zeta* (omega/ omega_n) ) A2 )  ; 

end 

%Identify  maximum  power  and  associated  design  point 
[x, y] =max (Power) 


A.6  Six- Variable  Model  Optimization  Code 
for  Factorial  Sampling 

%Optimal  Power  for  Six  Variable  Model  (Factorial  Sampling  Method) 
%Created  January  2015  by  Russell  Nelson 

%The  following  script  evaluates  the  average  power  dissipated  by  the  load 
%resistor  of  a  piezoelectric  generator.  The  equation,  consisting  of  6 
%independent  variables,  used  was  derived  by  Professor  Hong  Zhou. 
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clear  all 


a=l 

b=l; 

c=l; 

d=l; 

e=l  ; 

f=l; 

%Set  variable  limits (range) 

A_min_lim  =  0.01; 

A_max_lim  =  0.99; 

R_min_lim  =  .02; 

R_max_lim  =  100; 

Obase_min_lim  =  120*pi; 

Obase_max_lim  =  360*pi; 

C0_min_lim  =  0.00000001; 

C0_max_lim  =  0.000001; 

WSC_min_lim  =  120*pi; 

WSC_max_lim  =  360*pi; 
zeta_min_lim  =  0.005; 
zeta_max_lim  =  0.02; 

%Set  number  of  levels  per  variable 
levels=12 ; 

%Computes  values  of  power  for  all  combinations  of 
%variables  at  all  levels 

for  A  =  linspace (A_min_lim, A_max_lim, levels) 
b=l ; 

for  R  =  linspace (R_min_lim, R_max_lim, levels) 
c=l; 

for  Obase  =  linspace (Obase_min_lim, Obase_max_lim, levels) 
d=l; 

for  CO  =  linspace (C0_min_lim, C0_max_lim, levels ) 
e=l  ; 

for  WSC  =  linspace (WSC_min_lim, WSC_max_lim, levels) 
f=i; 

for  zeta  =  linspace ( zeta_min_lim, zeta_max_lim, levels) 


91 


P(a,b,c,d,e,f)  = (AA2*R*0base A  6* ( 1+CO  A2*0base A2  *RA2 ) A2 ) / . . . 
(2* ( ( (~l*ObaseA2+WSCA2) * (l+COA2*ObaseA2*RA2) + (C0*AA2) * . . . 
ObaseA2*RA2) A2+ (2*WSC*zeta*0base* ( 1+CO A2*ObaseA2*RA2 ) . . . 
+AA2 *Obase*R) A2)  )  ; 
f=f +1 ; 
end 
e=e+l ; 
end 
d=d+l ; 
end 
c=c+l ; 
end 
b=b+l ; 
end 
a=a+l ; 
end 

%Vectorizes  P  matrix 
Pvec=P ( : ) ; 

%Finds  the  maximum  power  in  Pvec 
[v, I ]  =  max (Pvec) 

%Determines  the  location  of  each  variable  level  for  maximum  power 
[j,k,l,m,n,o] =ind2sub (size (P)  ,  I)  ; 

%Determines  the  value  for  each  variable  at  maximum  power 
A  =  linspace (A_min_lim,  A_max_lim, levels )  ; 

Aopt  =  A ( j ) 

R  =  linspace (R_min_lim, R_max_lim, levels ) ; 

Ropt  =  R(k) 

Obase  =  linspace (Obase_min_lim, Obase_max_lim,  levels )  ; 

Obaseopt  =  Obase (1) 

CO  =  linspace (C0_min_lim, C0_max_lim, levels ) ; 

COopt  =  CO  (m) 

WSC  =  linspace (WSC_min_lim, WSC_max_lim, levels ) ; 

WSCopt  =  WSC (n) 

zeta  =  linspace (zeta_min_lim, zeta_max_lim, levels)  ; 
zetaopt  =  zeta(o) 
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A.7  Six-Variable  Model  Optimization  Code 
for  NOBLH  Sampling 

%6  Variable  Model  Optimization  Code  for  NOBLH  Sampling 
%Created  January  2015  by  Russell  Nelson 

%The  following  script  evaluates  the  average  power  dissipated  by  the  load 
%resistor  of  a  piezoelectric  generator.  The  equation,  consisting  of  6 
%independent  variables,  used  was  derived  by  Professor  Hong  Zhou. 

clear  all 

%Import  design  points 
filename  =  ' NOL6var . csv ' ; 

%Rename  file 
NOLH=csvread (filename) ; 

%Compute  power  for  all  design  points 
for  i  =1:512 

A=NOLH  (i, 1)  ; 

R=NOLH (i, 2) ; 

Obase=NOLH  (i, 3) ; 

CO=NOLH  (i, 4) ; 

WSC=NOLH (i, 5) ; 
zeta=NOLH  (i, 6) ; 

Power (1, i)  = (AA2 *R*Obase A 6* ( 1+C0 A2*ObaseA2*RA2 ) A2) /  .  .  . 

(2* ( ( (-l*ObaseA2+WSCA2) * (l+COA2*ObaseA2*RA2) + (C0*AA2) * . . . 
ObaseA2*RA2) A2+ (2*WSC*zeta*Obase* ( 1+C0 A2*ObaseA2*RA2 ) . . . 

+AA2 *Obase*R) A2) )  ; 

end 

%Identify  maximum  power  and  associated  design  point 
[ x , y ] =max (Power) 
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